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ABSTRACT 


/ 

The  need  for  tools  capable  of  handling  non-stationarities  in  the  spectral  content  of 

i'. 

the  data^ha?  been  recognized  as  early  as  1946.  The  Wigner-Ville  Distribution  (WD)  has 
been  extensively  used  since  its  introduction  in  1948,  but  suffers  from  some  associated 
problems  (e.g.,  spectral  cross-terms  and  requiring  the  use  of  analytic  signals).  An  alter¬ 
native  Distribution  is  proposed,  which  has  its  origin  in  the  definition  proposed  by  Page 
of  ^Instantaneous  Power  Spectrum"^  I  PS).  Its  characteristics  are  examined  and,  when 
pertinent,  compared  to  the  WD.  It  is  shown  to  be  less  sensitive  to  the  problems  aftieting 
the  WD.  but  provides  less  frequency  resolution.  The  usefulness  of  a  parametric  (AR) 
version  is  investigated.  Some  typical  test  signals  are  examined,  to  demonstrate  the  per¬ 
formance  and  trade-offs  of  IPS  and  its  parametric  version.  K-'-fi'  '  */•'  t  ^.  1  *  "  ' 
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I.  INTRODUCTION 


A.  THE  PROBLEM 

The  use  of  the  Fourier  Transform  as  a  spectral  description  of  signals  is  a  concept 
whose  usefulness  is  restricted  to  the  class  of  stationary  signals.  Though  mathematically 
elegant  and  convenient.  Fourier  decomposition  of  a  signal  can  often  mask  the  true 
spectrum,  since  the  assumed  basis-functions  implicitly  oppose  any  notion  of  time- 
dependency. 

When  the  spectral  content  of  the  signal  changes  with  time,  as  is  often  the  case  in 
fields  such  as  Communications.  Seismology  or  Speech  Processing,  more  powerful  tools 
are  needed.  The  time-dependency  of  the  spectral  content  should  be  apparent  and  meas¬ 
urable  through  the  use  of  a  more  general  type  of  representation  of  the  signal.  Gabor 
[Ref.  1]  proposed  such  a  representation  by  introducing  the  concept  of  a  spectral  de¬ 
scription  in  a  joint  time-lrequcncy  plane. 

Since  then,  several  attempts  have  been  made  [Refs.  2,  3,  4.  5]  to  derive  a.  or  the. 
function  capable  of  correctly  describing  the  distribution  of  the  signal's  energy  in  this 
plane.  Short  Time  Fourier  Analysis  (SI  FA)  has  been  widely  used  and  is  regarded  as  a 
very  convenient  approximation  to  the  “true"  distribution. 

However,  obtaining  this  distribution  is  an  under-determined  problem,  since  an  in¬ 
finity  of  valid  solutions  is  introduced  by  allowing  non-stationary  descriptions.  For 
example,  any  time  signal  can  be  represented  as  a  time-varying  DC  component.  To 
illustrate,  consider  the  following  discrete  time  signal: 


Two  equally  valid  representations,  in  the  sense  that  each  one  is  capable  of  generat¬ 
ing  the  observed  signal,  are  shown  in  Figures  1  and  2. 

Each  one  of  these  descriptions  explains  the  existence  of  a  sinusoidal  time-series. 
Though  there  is  no  conceptual  distinction  between  them,  the  algorithm  of  analysis  must 
have  the  ability  to  decide  which  one  will  be  taken  as  the  “true"  one,  based  on  an  ex¬ 
plicitly  or  implicitly  built-in  set  of  rules.  The  situation  thus  arises,  where  “one  will  see 
what  one  wants  to  sec". 


Figure  1.  x{n)  represented  as  a  constant  sinusoid 


Further  complications  appear  when  we  consider  the  fact  that,  given  a  time-series, 
we  cannot  determine  if  the  underlying  process  is  or  is  not  stationary,  without  using  ar¬ 
bitrary  and  often  inappropriate  assumptions  about  what  is  the  local  behavior.  When  the 
algorithm  of  analysis  produces  a  time-varying  spectrum,  the  fact  remains  that  there  is 
at  least  one  equally  valid  time-invariant  spectral  description  (the  Fourier  Transform 
based),  which  is  arbitrarily  ruled  out. 

B.  APPROACHES 

If  the  signal  is  to  be  represented  in  the  joint  time-frequency  plane  in  a  sensible 
manner,  the  distribution  must  have  some  key  properties.  For  example,  a  shift  in  the 
time-series  should  always  imply  a  corresponding  shift  of  the  spectral  representation 
along  the  time-axis.  Also,  a  multiplication  of  the  time-series  by  a  complex  exponential 


should  result  into  a  shift  of  the  spectral  representation  along  the  frequency  axis.  Without 
these  properties,  physical  interpretation  of  the  representation  can  be  an  impossible  goal. 
This  is  one  of  the  reasons  why  Ambiguity  functions  are  not  used  in  spectral  estimation. 

This  technique  of  imposing  constraints  that  are  felt  needed  in  a  distribution  can  be 
carried  further,  and  the  distribution  (or  one  of  the  distributions)  formed  in  this  refine¬ 
ment  process  defined  as  the  true  distribution.  One  of  the  most  successful  representations 
obtained  by  this  approach  is  the  Wigner-Ville  Distribution  (WD),  which  has  extensively 
been  used  since  its  introduction  in  1948  [Ref.  2]. 

A  different  approach  to  solve  the  indeterminacy  is  to  define,  a  priori  and  unambig¬ 
uously,  what  the  considered  true  distribution  is.  Though  this  approach  lacks  control  of 
the  properties  of  the  resulting  distribution,  it  has  the  advantage  of  more  closely  pre¬ 
serving  physical  meaning  in  whatever  t  ;sulu  it  produces.  A  typical  example  is  the 
definition  of  “instantaneous  Power  Soectrum"  proposed  by  Page  [Ref.  5],  uniquely 
determining  a  resulting  distribution  which  is  amenable  to  physical  interpretation. 

C.  OBJECTIVES 

The  Wigner-Ville  Distribution  has  been  extensively  studied  [Refs.  6.  7.  S.  9  ].  and  its 
characteristics  are  fairly  well  understood.  However,  iis  use  is  still  hampered  by  some 
problems.  First,  it  has  uneven  performance  for  different  classes  of  signals.  It  performs 
optimally  for  single-component  linear  FM  (linear  chirp),  but  has  worse  performance  for 
less  regular  spectral  dynamics.  Also,  when  multi-component  signals  are  present,  the  WD 
creates  artifacts  in  the  spectrum,  lying  mid-w«y  between  true  components.  These 
artifacts  arc.  up  to  some  degree,  recognizable  and  treatable  due  to  an  alternating;  sign 
[Ref.  6).  but  can  severely  mask  results  when  analyzing  more  complex  signals.  The  use 
of  analytic  signals  is  usually  required,  not  only  to  avoid  the  need  for  sampling  at  twice 
the  Nyquist  rate,  but  also  to  prevent  the  appearance  of  undesirable  artifacts  that  would 
otherwise  be  created  by  the  interference  between  positive  and  negative  frequencies 
[Refs.  6.  7.  10J. 

Though  it  has  not  received  proper  attention  in  the  literature,  Page's  Instantaneous 
Power  Spectrum  is.  once  properly  understood  and  treated,  a  practical  alternative  to  es¬ 
timate  the  time-varying  spectrum. 

The  study  and  development  of  such  an  alternative  is  the  main  goal  of  this  thesis. 
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II.  INSTANTANEOUS  POWER  SPECTRUM  (IPS) 

A.  DEFINITION 

In  an  attempt  to  accommodate  the  notion  of  instantaneous  frequency  content,  Page 
defined  the  Instantaneous  P^vver  Spectrum  as  the  derivative  of  a  running  energy  spec¬ 
trum  [Ref.  5]: 

P~(t.  J)  =  4^\S;(f\2  (2) 


where 


S7(J)  = 


(~'i 


That  is,  the  Instantaneous  Power  Spectrum  was  defined  as  being,  at  each  frequency, 
the  rate  of  change  of  the  energy  collected  by  a  Fourier  transform  taken  from  -  oo  .  up 
to  the  time  of  analysis.  This  concept  was  later  extended  by  Levin  [Ref.  3],  who  intro¬ 
duced  a  complementary  backward  run.  similarly  defined  ac: 


A*,  n  = 


•s,V) 


-i 

i 


Hi 


where 


roc 


A/)  = 


s(t)c 


J, 


and  defining  the  Instantaneous  Power  Spectrum  as  the  average  of  these  two  runs: 


iPS(,,J)  =  ±[p-(i.f)  +  p+(i.  /)]. 


(H 
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This  expression  can  be  put  in  the  form  [Refs.  3,  1 1  ] 

IPSO,  J)  =  Real  (7) 

where  S(/)  is  the  Fourier  Transform  of  s(r). 

IPS  can  thus  be  seen  to  be  the  real  part  of  Rihaczek's  distribution  [Ref.  4],  which 
allows  one  to  think  of  IPS(t.J)  as  being  the  energy  density  entering,  at  time  t,  an  infi¬ 
nitely  narrow  filter  centered  at  frequency/  An  enlightening  treatment  of  this  point  can 
be  found  in  Ackroyd  [Ref.  12]. 

I.  Time  Domain 

Let  us  now  consider  the  following:  we  want  a  function  of  /,  r  whose  Fourier 
Transform  is  IPS(i.  f)  .  That  is,  we  want  G(r.  t)  such  that 


IPS(Lf) 


G{[,  i)e~^"^dr  . 
J 

— :>c 


(8) 


Then,  using  (7) 


Gil.  T) 


IsiOSXfte-^'  +  s(/)*S(/)e/2*/']  e^V/ 


(9) 


from  where,  with  the  appropriate  variable  substitutions  we  obtain 


G[i.  r)  =  ~ 


/» 

r- 

S{! ) 

5  (/  -  tj) e~i“'J,1dtj  +  s'(i) 

s(i  +  }j)c~j2~fnd)] 

- 

JW 


df.  (10) 


Interchanging  the  order  of  integration,  we  get 


GO.  ')  =  4*  -  t)  +  5*(/)s(r  +  t)]  . 


(11) 
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IPS  can  now  be  res’ated  as: 


OO 

r 


IPS{i,f  =  ~ 


[s(/)s  {t  —  r)  +  s*(t)s{t  +  t)]  e  J2~J'd~  . 


—CO 


(12) 


2.  Frequency  domain 

Alternatively,  a  dual  expression  can  be  derived,  expressing  IPS  in  terms  of  the 
Fourier  Transform  of  the  signal. 


IPS(t,j)  =  -j 


s(t)s  (t  —  r)e  j2"^dr  +  ~ 


s  ( i)s(t  +  r)e  "J'dr 


J 

— OO 


oc  oc 

r 


s^^'dn 


e~j-f'-dx  + 


—  OC  — OC 

OC  oc 


—  OO 


*  r 

S  (u) e~j2"“rdfu 

J  V 

S(y)ca,:y<J+r)dy  e'J2'frdx 


—  OC  -X 

OO  cc- 


Sl^sUyy^e-^'Siy  -j)dydu  + 


%) 

— OC  — oc 


(13) 


+  T 


5  (n)S(y)e  j2rix‘  ei2"y,d{y  -fjdyd^i 


J  J 

—  OO  — OO 


S{n)S{f)e-J2r:(f-rt'dn  +  -i- 


J 

— OO 


S‘(M)S(/)c-y2fr(//"/)W  • 


— OO 
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With  the  appropriate  variable  substitutions,  we  get 


IPSb,  J)  =  T 


S(f+  y)S*(f)ell~y'dy  +  ± 


S'(f-y)S(f)e,l7:vdy  (14) 


%> 

— OC 


and,  finally 


oo 

r 


irs{t.j)  =  y 


[S(/)SV-  V)  +  S*{f)S[f+  y)]  en*1,dy  . 


— OO 


(15) 


3.  IPS  for  discrete  signals 

The  discrete  version  of  IPS  follows  directly  from  (  12  )  as 


IP  Sin.  (?)  = 


AT 


zLi 

A'=—  >c 


[s(/;)s  («  -  k)  +  s*(«)s(«  +  Clje  ^/\ 


(16) 


A  formal  derivation  of  (  16  )  can  be  found  in  Appendix  A. 

B.  RELATIONS  W  ITH  OTHER  TIME-FREQUENCY  DISTRIBUTIONS 
1.  Rihaczek's  Distribution 

The  Rihac/ek  Distribution  was  proposed  in  196S  as  the  true  energy1  distribution 
in  an  attempt  to  unify  existing  results. 

Formally,  the  distribution  is  [Ref.  4] 

t{t.  J)  =  s(i)S'(f\e~J2'f‘  (17) 


where  s(t)  is  the  analytic  signal  being  analyzed,  and  S(/)  is  its  Fourier  Transform. 

As  was  derived.  dr0.  f„)  is  the  complex  energy  density  at  the  point  in  the 

time-frequency  plane.  A  similar  expression  had  been  considered  by  Levin,  who  defined 
the  Complex  Instantaneous  Power  Spectrum  as  the  complex  conjugate  of  (  17  ). 


Rihaczck  used  the  concept  of  "complex  energy",  whose  real  part  is  the  real  energy,  while  the 
imaginary  part  is  the  reactive  energy. 
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Using  Levin's  notation,  the  Complex  Instantaneous  Power  Spectrum  is 


AXt.J)  =  s\t)S{f)e,2*fi,  (IS) 

where  s(t)  is  the  complex  envelope  of  the  real  signal  under  analysis. 

We  thus  see  from  (  7  )  that,  though  addressing  a  broader  class  of  signals,  IPS 
is  both  the  real  part  of  Rihaczek's  Distribution  and  the  real  part  of  Levin's  Complex 
Instantaneous  Power  Spectrum.  Formally, 

IPS(t,J)  =  Real[c(/,  /)]  .  (19) 

2.  Cohen's  Generalized  Phase-Space  Distribution  Functions 

In  1966.  Cohen  introduced  a  generalized  class  of  time-frequency  signal  repres¬ 
entations,  given  by 


oc  oc  oc 


c{r.  J)  = 


O(o,  r)  s(tx  +  -^-)/(/,  -\)  e12'*"'  Ll~'J)dvdt,dr  (20) 


J  J  J 

— DO  — OO  — OO 


where  the  choice  of  cf>(i),  r)  will  determine  the  resulting  distribution  |Ref.  l?j. 
If  we  slightly  rearrange  (  20  )  to  read 


OC  >7  OO 


(% 

- 

/% 

C(t,  J)  = 

%> 

0(u.  t) 

— *J 

s[t\  +  y  )s*(r,  -  *7  )c,2'l!ldtx  e 

— OO  — OC  — OC 


t+-/) 


dod T  , 


(21) 


we  can  interpret  it  as  the  double  Fourier  Transform  of  T(i>,  r)  [Ref.  14],  where 


OO 


i  • 


J 

—so 


(22) 


We  can  now  recognize  (  22  )  as  a  general  expression  for  the  well  known  Ambi¬ 
guity  function,  which  can  serve  as  a  genera!  definition  of  the  time  frequency 
autocorrelation  function  [Refs.  3.  4  ). 


Each  choice  of  weighting  function  <£>(t>,  t)  in  (  20  )  provides  a  different  distrib¬ 
ution,  because  each  one  defines  the  combined  autocorrelation  (  22  )  in  a  somewhat  dif¬ 
ferent  way. 

Let  us  see  which  distributions  result  from  three  different  realizations  of  r) 
(complex  exponential,  constant,  and  real  sinusoid)  [Refs.  4,  13,  14  j. 


a.  Complex  exponential 


Ofu,  t)  =  e  2 


The  combined  autocorrelation  becomes 


(23) 


T(i-.  r)  =  c'"1 


or- 

r 


*0i  +  y  )*’(',  -  y  ) 


— OC 


r 

Am  Am 

— oc 


and.  defining  a  new  variable  /  = 


/,  +  y.  we  obtain 


(24j 


oc 

/* 

'Efe.  t)=  5(/,)s*(f,  -  .  (25) 

* 

—  OC 


The  resulting  distribution  is 


OC  OC  OC 


J)  = 


s(r,)s  (r>  —  r)e,'lnU>  e  J2n^I+'^d[,dTdu  , 


(26) 


—oo  — oo  — oc 


Q 


C{t,  f)=  \  1  s{t\)s\h  -  r)e  j2'r/d(t1  -  t)dt^dx 


— OO  — OO 

oo 


(27) 


=  |  s(t)s  ( t  —  r)e 


=  5(/)S  {J)e 


which  is  the  well-known  Rihaczek  Distribution. 

Though  this  definition  of  the  combined  autocorrelation  is  probably  the 
most  intuitive  one.  its  lack  of  mathematical  symmetry  has  the  undesirable  effect  of 
making  the  distribution  complex.  It  has  extensively  been  used  in  radar  theory,  since 
(  25  )  is  exactly  Woodward's  definition  of  the  Ambiguity  function  [Ref.  3]. 


b.  Constant 


<h(t>,  t)  =  1 

I  The  combined  autocorrelation  becomes 


(28) 


mr)  = 


s(h  +  Mb  -  *T  )  el2~l,'dt) 


4/ 

—  DC 


(29) 


The  resulting  distribution  is 


OO  oc  oo 


c(t,  J)  = 


s(t]  +  ~)s  (f,  — ^-)c/2  ‘lf‘  e  ^"lu+'^dt^dTdo  , 


(30) 


— OO  — oo  — oo 
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OO  DO 

r  r 


cU,  J)  = 


s(f  +  rf-  )s\r{  -~r)e  j2r-'-fd{ix  -  t)di xdt 


— OC  — OO 
OO 

r 


(31) 


s(t  +  J  )s\t  -  Y  )e  J2*VdT 


— oe 


which  is  the  Wigner-Ville  Distribution.  This  definition  of  the  combined  autocorrelation, 
though  only  slightly  different  from  the  one  in  (  25  ),  has  the  needed  symmetry  to  make 
this  distribution  real,  and  mathematically  very  convenient  (when  dealing  with  analytic 
signals). 


c.  Cosine  function 


<l)(o,  t)  =  cos(ttdt) 


The  combined  autocorrelation  becomes 


(32) 


T(o,  t)  =  cos(-I)t) 


•S(fi  +  -3j-  )5  (/,  -  -^  )t^ri  V/, 


+  e~J™) 


•s(fi  +  )5*(fi  — 7"  )ci2~l!,dt] 


%/ 

— OO 


(33) 


— oc 


+ 


s(t]  +  y  )s\t]  —  ~  )e,2~'j(''  ~  "i"  )dtx 


J 

— oe 
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we  get 


Substituting  t  for  in  the  first  integral,  and  i  for  q  —  in  the  second, 


oc  oo 


r- 

f* 

_ 

t(M)=j 

s(r)s  (i  —  r )e*2™'dt  + 

5  (/).s(/  +  d)/1™1  dt 

—4/ 

* 

-oo  ~ 

-oc 

oo 

I* 

~  [s(/)s*(f  ~  T)  4-  s*(/)5('  +  T^^'dl . 

At 

* 

—  OO 


(34) 


The  resulting  distribution  is 


OO  OO  OC 


r  r 


c(f.  J)  = 


~  -  r)  +  s\r])s(tl  +  T)y™'e 


2m,-j2~(ti+zJ) 


dixdvdr 


-oo  — OC  —  OO 
OO  OC 


4"  -  t)  +  s*(/i)s(/,  +  t)]e  J2':fd(il  -  (35' 


— oo  —oo 


[s(7)s  (/  —  t)  4-  s'(l)s(l  +  t)]c  j2~'Jdr 


which  is  the  expression  we  had  found  earlier  for  IPS,  In  this  case,  the  definition  of  the 
combined  autocorrelation  also  possesses  the  symmetry  needed  for  a  real  distribution. 
The  kernel  that  generates  IPS  (  0(l>,  t)  =  cos(7u>t)  )  has  been  considered  by  Cohen 
(Ref.  13J.  and  shown  to  generate  the  Margeneau-Hill  distribution  [Ref.  15).  well  known 
in  quantum  mechanics  theory. 

3.  Wigner-Ville  Distribution 

Since  the  Wigner-Ville  Distribution  (WD)  has  enjoyed  wide  acceptance,  we  will 
define  more  closely  the  relations  between  IPS  and  WD.  As  an  indirect  connection,  we 
can  relate  both  to  the  Rihaczek  Distribution.  From  (  29  )  the  WD  implicitly  defines  the 
combined  autocorrelation  as  [Ref.  4] 


oc 


Y(of  r)  = 


S(t]+^)s\il-^)e/2ru'c1i1 


=  e~Jr'^ 


s(t:)s\i i  -  r)e>'lrilt]dtl  . 


Hence,  the  Wigner-Ville  Distribution  is  given  by  [Ref.  4] 


OC 


- 

r  -i 

H  D(t,  f  =  Ft  . 

e-J™ 

L-  % 

i 

V 

t- 

1 

* 

—  DO 


(36) 


(37) 


where  FL,.[  ■  ]  denotes  the  2-D  Fourier  transform.  From  the  convolution  property  of  this 
operator,  and  using  (  25  ), 


IF  DU.  J)  =  Fv **  rt 


OC 

r 

j(/,)^*(/i  -  r)e,'2ri‘ld!l 

—  DC 


=  FvJj>'J™]  **c(/.  f] 
=  eiirfl  **  di,  J) 


(38) 


where  **  stands  for  2-D  convolution,  and  c(i,  f)  is  the  Rihaczek  Distribution.  From 
(  3S  ).  using  (  Id  )  and  the  realness  property  of  WD  (see  Section  II-C).  we  obtain 


lVD(i.  f)  =  Real[  e>A~fl  **  e(i.  /)] 
IPS(t.  f)  =  Real[  6(i,  J)  **  e(r,  /)]  . 


(39) 


We  thus  see  how  the  two  distributions  can  be  obtained  from  the  Rihaczek  Distribution 
by  an  appropriate  choice  of  the  function  convolved  with  c(t,  J) . 
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We  found  earlier  that  each  choice  of  ^(o,  r)  in  Cohen's  class  of  Distributions 
(  20  )  yields  a  particular  definition  of  the  combined  autocorrelation  function 


H>,T)  =  0(D,  r) 


%} 

— OO 


(40) 


and  that  the  2-D  Fourier  transform  of  (  40  )  resulted  in  a  particular  distribution,  de¬ 
pending  on  O(o,  t). 

We  see  now  that  (  40  )  is  the  product  of  two  functions,  $(i>,  r)  and  the  integral. 
Its  double  Fourier  Transform  is  thus  a  double  convolution  in  the  transform 
domain. [Ref.  14] 

That  is: 


c(t,j)  =  FVixnv,  t)] 


o© 


K,-. 

0(0,  T) 

s{t\  +  "T  )s*Ui  —  \  )ej7~u'dt] 

—  %> 

I  ** 

=  y,[0(u,r)]**ftiT 

*(', +y)*V, -t)^V;, 

X 

=  Ft  .[<I>(o.  Tl]  **  IVD(t.f). 


In  the  particular  case  of  IPS, 


JPS{i ,  f)  —  Fv  .[  cos(7tut)]  **  \YD{i. f) . 


(42) 


Therefore. 


IPS{t.  f)  =  cos(4r//)  **  U  D{t,  f) ,  (43) 

which  we  could  have  deduced  directly  from  (  39  )  and  the  fact  that  the  \YD  is  always 
real,  as  follows: 


IPSO,  J)  =  Real[c(f.  /)] 

=  y  W'.  J)  +  t\i,  /)] . 


(44) 
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Using  (  38  ), 


IPS(i ,  f)  =  y  [e-74^'  **  IFZ>(/,  y)  +  ^4rr/l  **  lVD(i.  /)] 

=  y  [e-7^'  +  e7'4*7']  **  HZ)(/,  /) 

=  cos(4?r//)  **  lFZ)(f,  y)  . 


(45) 


At  this  point,  we  would  like  to  be  able  to  express  the  WD  in  terms  ofIPS,  with 
a  relation  of  the  form: 


\VD{t,  f]  =  OUPSU,  rn , 


(46) 


where  0[ .  ]  is  an  arbitrary  operator.  Unfortunately,  the  inverse  filter  does  not  exist, 
since 


1 

COS(  7TUT ) 


(47) 


is  not  defined  over  the  required  range  of  n  and  r  . 

4.  Short-Time  Fourier  Transform 

Though  the  Short-Time  Fourier  Transform  (STFT)  can  not  be  considered  a  true 
Time-Frequency  Distribution  [Ref.  14],  its  widespread  use  fully  justifies  the  study  of  how 
it  relates  to  IPS.  This  hs  been  addressed  in  [Ref.  11],  and  will  be  done  here,  after  which 
some  indirect  results  will  emerge.  Uet  us  thus  consider  the  STFT  as  defined  in  (  48  ): 


STFT{t,  to)  =  |F(t,  co)  | 2 


s(/  +  t)iv  (t)( d^'dr 


J 

—  OC 


s(t  +  r)s  (/  +  o)vv  (t)w(i))c  e1™  dzdv 


(4S) 
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Substituting  j(r  +  r)  and  vv"(r)  by  their  Fourier  definitions,  we  have  that 


GO  OO 


STFT(t,  to)  = 


s  (t  +  o)vv(u) 


S{^(,+z)dz 


— OC  — 'OO 
OO  OO 

r  r 


— OO 


\V  (p)e  'Jp'dpe  dxdv 


%/ 

— OO 


s  (/  +  u)vv(u)S(p  +  io)lV  {p)^u3+p^'e^u>vdvdp 


(49) 


— oo  —oo 
OO  oo 


s  (r  +  o)S(p  +  co)^(p+“j>(,+l>w(u)H/  (p)e  ^pvdvdp  . 


— OC  — oo 


We  thus  see  that  the  STFT  is  a  time-frequency  smoothed  version  of  the  Rihaczek  dis¬ 
tribution  of  the  signal.  The  smoothing  function  is  the  Rihaczek  distribution  of  the 
window  used  to  compute  the  STFT.  Since  the  right-hand  side  of  (  49  )  is  always  real  and 
positive,  and  the  previous  relations  are  valid  for  arbitrary  signals  and  windows,  we  con¬ 
clude  that  windowing  a  Rihaczek  distribution  with  a  Rihaczek  distribution  guarantees  a 
real  and  everywhere  positive  distribution.  This  result  will  be  important  when  addressing 
the  issue  of  positivity. 

C.  BASIC  PROPERTIES  OF  IPS 

As  was  already  mentioned,  and  is  easily  seen  from  (  20  ).  it  is  the  choice  of  a  par¬ 
ticular  O(o.  t)  that  determines  the  particular  distribution  and.  hence,  its  properties.  It  is 
thus  desirable  to  establish  direct  connections  between  the  properties  of  Olr.  r)  and  the 
properties  of  the  resulting  distribution.  This  issue  has  been  addressed  |Ref.  14).  and  we 
will  in  the  most  part  only  state  results,  as  applicable  to  IPS.  The  basic  properties  of  IPS 
are: 

1 .  Time  shift 

If  aft)  =  s(t  -  /„).  then 

IPSJt.J)  =  IPSs(,-t0.f).  (50) 

2.  Modulation 

If  aft)  =  s(t)e2‘ o',  then 

IPSjt,f)  =  IPS5(t.f-fv).  (51) 
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3.  Marginal  in  Time 


OC 


IPS(t,  f)  df=  |  *(/)  i2  . 


—  OC 


4.  Marginal  in  Frequency 


OC 

lPS{t .  f)dt=\  S\J)  |2 . 
J 

— OC 


where  S(/)  =  F[s(r)]  . 

5.  Realness 

lPS(t.  J)  is  real  for  all  /,  /. 

6.  Instantaneous  Frequency ‘ 

OC 

/•/«(/.  ./)# 

c 

— - t - =  yj( l)  (instantaneous  frequence) 

I  .'( / )  i' 


7.  Group  Delay 

S 

r 

t  iPSu.  'f)di 

— - - ; - =  TJ  f)  (group  delay) 

i-Vir 

S.  Zero  Power 

sU.)  =  0  =>  IPS(r0.  J)  =  0 
5(/,|  =  o  =*->  /c)  =  0 


(52) 


(53) 


(54) 


(561 


(5") 


The  Wigncr-Villc  Distribution  also  respects  properties  1-7  [Ref.  6],  but  only  a 
weaker  version  of  property  8. 

D.  FURTHER  PROPERTIES 

Though  the  properties  presented  in  Section  II-C  constitute  the  basic  set  of  rules  to 
which  IPS  is  bound,  some  further  anahsis  must  be  made.  The  effects  of  linear  operators 
should  be  investigated,  to  fully  understand  the  behavior  of  IPS. 


This  property  is  valid  for  anahtic  signals  onh  (  see  .Appendix  B  ). 
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1.  Windowing  in  the  time-domain 

The  issue  here  is  to  determine  how  IPS  is  affected  by  windowing  the  original 


time  signal. 


We  will  at  this  point  introduce  the  notation: 


IPSli,  J)  =  j-  IPSj{t,  y)  +  y  IPSpt,  f) 


where 


IPSd  U>  f)  =  d(t)d  (r  -  t)c  u  '  'dx 


IPSpt ,  J)  =  /(/)</(/  +  x)c~J’2rP'dx 


If  the  signal  is  windowed,  that  is: 


d(i)  =  s(/)u  (/) . 


then  the  IPS  of  the  windowed  sienal  will  be 


IPSjUf)  =  P  [d,U)d  (t-x)+d  (t)d(i  -  r)]a_y*"Vr 


=  —  [s(/)vv(/)s  (i  -  r)iv  (/  -  r)  + 


4-  5  (/)»■  (r)s(r  +  t)w(i  +  t)]c  dx  . 


Separating  the  two  terms. 


IPSaU,  /)  =  -7 


[s(r)i*(f  -  t)][w(;)iv*(/  -  r)>  J2"f'dr 


—  X) 


+  ' 


[5  ( t)s(i  +  T)][w  (r)w(r  +  t)]c  j2"^’dr 


— OC 
CkC 


F  \lPS~]F  \lPSZ\e~nKfxdx 


F  1  UPS*]  F  \lPS*.]e~jlrJ''dT 


IPS s  (t.  v)IPS w  (/,/-  v)dv  4-  — 


(63) 


IPS*(t.  u)IPS*(i,f-  v)dv  . 


J 

~3C 


where  F  [.]  is  the  inverse  Fourier  Transform  operator.  Hence 

IPSJI.  /)  =  4-  [ips;u.  J)  *  IPS;.{t .  /)]  +  y  [lPS*(L  f)  *  7/>S+(r.  f) 


(64) 


where  *  denotes  eons  olui ion  in  the  f  variable. 

That  is.  while  1'or  the  Wigncr-Ville.  windowing  the  time  signal  implies  the  con¬ 
volution  along  the  frequency  axis  of  the  signal's  WD  with  the  window's  \YD  (Ref.  6j. 
for  IPS  we  have  a  sum  of  two  convolutions:  the  first  term  of  the  signal's  IPS  (IPS;)  is 
convolved  with  the  first  term  of  the  window's  IPS  (IPS;),  and  a  similar  convolution  is 
applied  to  ( IPS: )  and  ( IPS;)  . 

A  special  case  occurs  when  the  window,  the  signal,  or  both,  possess  conjugate 
symmetry  about  t,  the  time  at  which  we  are  computing  IPS  (  i.e., 
s(i  trl  =  s'(  1  —  t).  wit  +  r)  =  w'(t  —  t),  or  both  ).  If  w(f)  possesses  that  symmetry,  then 

IPSjt,  J)  =  IPS^d,  j)  =  IPSJi,  J)  (65) 


and,  by  linearity  of  the  convolution  operator.  (  64  )  becomes 


OO 

/» 


IPSJb,  f)  = 


lPSs(t,  v)IPSw(i,  f-  u)d V  . 


— OO 


(66) 


We  thus  see  that  in  this  case  IPS  and  WD  respond  similarly  to  the  windowing 
of  the  data:  the  IPS  of  the  windowed  signal  is  the  convolution  (in  frequency)  of  the  IPS 
of  the  signal  with  the  IPS  of  the  window. 

In  summary  ,  the  IPS  of  a  windowed  signal  is  a  frequency-smoothed  version  of 
the  IPS  of  the  unwindowed  signal.  The  time-resolution  is  not  affected. 

2.  Convolution  in  the  Time-domain 

The  question  now  is:  how  is  the  IPS  of  a  signal  affected  if  the  signal  is  pre- 
processed  with  a  filtering  operation?  So.  let 


d{i)  =  s(0  *  h(t) . 


(67) 


Then. 


IPSJt.  f)  =  Real  [OV)  *  /i*(/)]  TOO)  *  /t(/)]  eiW‘] 

=  Real  [is\t)  *  h'(t)']S{J)H[/)c/2 ~-;7] 

=  Real  [[/(/iSt/V'2"'7]  *  [/i*(t)//(/)r;2r/‘]]  , 


where  we  used  the  fact  that 


or 

/* 


e/2'/7[a(0  *  MO]  =  ^ 


a(r)b{l  —  x)dr 


'-DO 


—  r)e^2"^‘  ‘Wr 

=  *  [MOr^7']  . 


(69) 


We  thus  see  that  (contrary  to  what  happens  in  the  WD3  )  filtering  the  signal 
does  not.  in  general,  correspond  to  the  convolution  in  time  of  the  signal's  IPS  with  the 
IPS  of  the  impulse  response  of  the  filter.  What  is  observable  from  (  68  )  is  that  the  IPS 
of  the  filtered  signal  is  now  the  real  part  of  the  convolution,  along  the  time  axis,  of  the 
two  Rihaczek  Distributions.  The  net  effect  pi  educed  on  IPS  is,  thus,  still  a  smoothing 
operation  along  the  time  axis.  The  frequency  resolution  is  not  affected. 

3.  IPS  and  Moral's  Formula 

As  is  often  stated  in  the  literature,  IPS  does  not,  in  general,  respect  M oval's 
formula  [Ref.  16].  When  applied  to  the  WD,  M oval's  formula  states  [Ref.  17]: 


WDs{t.  J)\VDg{i,  j)didf=  <  s,g  >  <  s.g  > 


where  <  s,g  >  is  the  inner  product  of  s  and  g  .  Let  us  now  see  what  happens  in  the  IPS 


ipss it.  fjirsji.  f)dnif= 


rjr  or  or 


[s(rb  (i  -  t)  +  s  ( i)s(r  +  x )]e  jl'!''dx 


[£(/’)£*(/  -  o)  +  g\l)g(t  +  u)]f  j2~fldv  didf 


The  WD  of  a  filtered  signal  is  the  convolution  along  the  time-axis  of  the  original  signal's 
WD  with  the  WD  of  the  impulse  response  of  the  filter  [Ref.  6). 
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oo  oc 

r  r 


ips5(i ,  nipsj/,  f)didf= 


— oo  — oc 


oo  oo  oo 


0(/)5  (/  -  t)  4-  s  (l)s{t  +  t)]1>(/)£  (/  -  o)  +  g  (l)g{l  +  d)] 


—oo  — oc  — oo 


-flcAi+i) 


drdtdo 


(72) 


oo  oo 


4 


MO*  (/  -  r)  +  5  (/)*(/  +  r)][>(0g  (/  +  r)  +  g  ( t)g(t  -  r)]  di  dr 


— oc  — OC 


=  —■  <s,g>"  <s.g>  +  <s,g>  <s.g> 


+  ~  Real 


oc 

oc 

/* 

s(t)g{r ) 

s*(t  -  r)g\l  +  r)dr 

dt 

% 

-oo 

-oc 

and,  finallv 


OO  'OC. 

r  r 


/PS ,(t.  f  )IPSJj,  f)didf= 


%J  %> 

— oc  — oc 


-7  <s,  g>  <  s,g >  *  4-  Real 

L 

oc 

OC' 

{+ 

s(r  )g(0 

'  s\t  -  r)g\t  4-  T)dr 

dl 

-oc 

V 

OC 

(73) 


The  result  confirms  that  M oval's  formula  is  not,  in  general,  valid  for  IPS.  Only 
in  the  trivial  case  of  one  of  the  functions  being  a  real  constant  would  the  second  term 
in  the  right  hand-side  of{  73  )  reduce  to  -7  <s.g>  <s,g>\  and  hence  would  IPS  obey 
M oval's  formula. 
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4.  Recovery  of  time-signals 

To  address  he  problem  of  recovering  a  signal  from  its  IPS.  we  will  find  it  con¬ 
venient  to  consider  two  separate  cases. 

a.  Infinite  duration  signals 

For  reai  signals  with  5(0)  ^  0,  recovery  of  the  signal  is  always  possible  ex¬ 
cept  for  a  sign  indeterminacy,  by  considering 

IPS{t,  0)  =  Real[s(0S(0)] 

(74) 

=  S(0W) 

and,  since  |  5(0)  |  can  be  found  by  using  Property  4  (  53  ),  only  the  sign  of  5(0)  remains 
undetermined.  No  similar  results  are  available  for  complex  signals  of  infinite  duration. 
However,  if  s(r)  possesses  conjugate  symmetry  around  some  r0.  then  we  can  also  recover 
complex  signals  of  infinite  duration  up  to  a  phase  term,  by  considering  F~]  [  IPS(ta.J)  ]. 
Similarly,  if  S(f)  possesses  conjugate  symmetry  around  some  f  and  5(0)  #  0,  we  can 
again  recover  complex  signals  of  infinite  duration  up  to  a  phase  term,  by  considering 
F[  lPSU,f0)  ]. 

b.  Finite  duration  signals 

For  finite  duration  signals,  recovery  of  real  or  complex  signals  is  always 
possible  except  for  a  constant  phase  term  (or  a  sign,  in  the  case  of  real  signals). 

Let  us  denote  the  starting  time  by  t  =  0.  From  (  12  ) 

F~\lPS( 0.  /)]  =  \  [  s(0)s'(  -  t)  +  5*(0)s(t)  ]  (75) 


and,  since  s(  —  t)  =  0  for  r  >  0  . 

F~\lPS{0.  m  l->0  =  y  5*(0)s(t)  (76) 

which,  by  property  3  (  52  ),  leases  an  unknown  phase  term. 

E.  IPS  -  A  GENERALIZATION  OF  THE  WIENER-KHINCHIN  THEOREM 

By  redefining  auto-correlation.  IPS  offers  a  very  straightforward  extension  of  the 
Wiener-Khinchin  theorem  to  the  non-stationary  case,  as  we  shall  see. 

Let  us  define  the  auto-correlation  function  of  a  signal  as 

R(t.  r)  =  -y  E{s(t)s  (t  -  r)  4-  s*(i)s(t  +  t)]  .  (77) 
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We  have,  then,  using  (  12  )  and  by  linearity  of  the  convolution  operator  the  follow¬ 
ing  pair  of  relations: 


oc 

r 


E{fps(i,  m  = 


-oo 


R(t.  x)  e-^dx 


oo 

r 


R(t,  t)  = 


— oo 


E{IPS(t,  J)}  e**df 


These  relations  can  be  read  as  a  generalization  of  the  Wiener-Khinchin  theorem  to 
non-stationary  signals,  since  for  stationary  signals  they  promptly  reduce  to  that  theorem. 
If  s(t)  is  a  stationary  signal,  then 


R{r,  t)  =  y  £[s(/)s*(/  -  r)  +  s\t)s{t  +  t)} 

=  y  £{s(/)s‘(/  -  t)}  +  ^  £{$*(*)*(*  +  t)J 
=  Y^(t)  +  Y/?‘(-T^ 

=  £(t) 


E{IPS(t.J)) 


OC 

f  R(x)e~J2~S~ dx  =  |  S{J)  |2 


— OC 


(79) 


The  expected  value  of  IPS  is,  thus,  for  stationary  signals,  the  usual  Power  Spectral 
Density  (  PSD  ). 

Since  in  the  non-stationary  case  we  cannot  try  to  infer  ensemble  averages  from  time 
averages,  the  IPS  of  a  signal  is  one  realization  of  the  generalized  Wiener-Khinchin  re¬ 
lations  of  the  underlying  stochastic  process. 

A  different  view  of  IPS  when  applied  to  stochastic  processes  can  be  found  in  Grace 
[Ref.  11].  where  E{lPS(t,  J)}  is  shown  to  be  the  Fourier  Transform  of  Loeve's  General¬ 
ized  Power  Spectral  Density. 
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III.  RETHINKING  IPS 


A.  ON  THE  USE  OF  WINDOWS 

IPS,  as  given  in  (  12  ),  suffers  from  intense  ringing  at  all  frequencies  “touched”  by 
the  signal.  This  same  effect  is  one  of  the  reasons  why  Rihaczek's  Distribution  did  not 
find  much  acceptance.  An  example  can  be  seen  in  Figure  3,  where  a  linear  chirp  was 
used  as  the  test  signal. 


Figure  3.  IPS  for  a  linear  chirp 


This  ringing  effect  can  easily  be  understood  by  slightly  rearranging  (  12  )  to  : 
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With  the  appropriate  changes  of  variable,  we  get 


lPS(t,  f)  = 


1 


s(/) 


5\h)e-JW~  h)du  +  s*(0 


— oo 


=  y  |>to  WO  *  ^Ix1  +  s\t)  [s(0  *  e/2c/']] 
=  Real  [  5*(0  [5(0  *  ]  • 


(81) 


But  5(0  *  e>7’F'  is  the  output  of  a  filter  whose  impulse  response  is  a  complex  exponential. 
IPS  is,  thus,  for  each  frequency,  the  real  part  of  the  product  of  the  signal's  complex 
conjugate  and  the  output  of  an  infinitely  narrow  non-causal  filter  (non-causa!  oscillator) 
(Figure  4).  [Ref.  12] 


Figure  4.  Model  of  IPS 


The  ringing  is  hence  unavoidable,  and  will  persist  until  the  filter's  impulse  response  be¬ 
comes  negligible.  We  should  consider  stable  filters,  with  decaying  impulse  responses,  if 
we  want  to  diminish  this  effect.  This  reasoning  leaves  us  very  close  to  the  work  of  Fano 
[Ref  18],  which  was  later  extended  by  Schroeder  and  Atal  [Ref.  19]. 
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If  we  replace  the  oscillators  in  the  IPS  definition  by  stable  filters  with  impulse  re¬ 
sponses  of  the  form 

h(t)  =  ,  (S2) 


where  w(r)  is  any  real  and  symmetric  function  of  r  with  general  lowpass  characteristics 
and  w(0)  =  1,  then  the  resulting  distribution  D(t,  f)  is  given  by 


Ds(t,  J)  =  Real  [s'(t)ls(t)  * 

=  \  +  s\t)[s(t)  *  w(/)^r//]] 


(S3) 


=  4*  s(/)s  (r)u-(t  -  x)e  ;2rr^'  )dx  +  ~ 


s'(t)s{u)w{t  - 


Substituting  /,  for  (r  —  t)  in  the  first  integral,  and  — for  (t  —  0)  in  the  second  integral, 
we  get 


DSU.J)  =  t 


[s(t)s'(t  -  /,)  +  s’(r)i(/  +  /,)]  u-(/,)  e 


(84) 


where  y(r)  =  iv(t  — r). 

That  is.  the  resulting  distribution  is  the  IPS  of  the  windowed  signal,  where  the  win¬ 
dow  is  the  envelope  of  the  filter  impulse  response,  centered  at  the  time  of  analysis.  Since 
u(/  -  t)  is  real  and  symmetric  around  1.  we  see  from  (  66  )  that  a  smoothing  of  the  IPS 
in  the  frequency  direction  will  result. 

IPS  in  general  can  hence  be  written  as 


OO 

IPS(t,  /)  =  y  J  [s(/)s  (i-t)  +  s  { t)s(t  +  t)]  vv(0)u-(t)c_-/2'/'<Yt  .  (85) 

—OO 


where  u(r)  is  the  chosen  window. 


To  illustrate  the  effectiveness  of  the  windowing  operation,  we  computed  the  IPS  of 
the  test  signal  used  in  Figure  3.  A  Hamming  window  was  used  in  (  85  ),  and  the  results 
are  shown  in  Figure  5.  As  is  observable,  the  ringing  is  avoided. 
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Figure  5.  IPS  for  a  linear  chirp  -  Hamming  window 

B.  PROPERTIES  AFTER  WINDOWING 

Some  of  the  properties  of  IPS  will  be  affected  by  the  windowing  operation.  It  will 
thus  be  necessary  to  determine  which  properties  are  modified  and  establish  the  con¬ 
nections  between  the  characteristics  of  the  used  window,  and  the  effects  it  generates. 

Tracing  back  the  windowing  operation  to  the  kernel  function  (Ofa,  r))  in  (  20  )  that 
it  implies,  we  will  be  able  to  use  well  established  results  [Ref.  14]  about  these  functions. 
The  resulting  kernel  function  will  be 

0(i>,  t)  =  w(0)vv(t)  cos(7u>t)  ,  (86) 

which  can  be  proved  using  (  19  )  as  follows  : 


oo  oo  oo 


n>s(t,j)  = 


w(0)w(t)  cos(7tut)s(/j  4-  —  )s  (/j  — --  )(^lnvt'e  Jln<'v,+'^dvdrdtl  .(87) 


— oo  — oo  — oo 


2S 


Replacing  the  cosine  by  complex  exponentials,  we  have 


oo  oc 


IPSU,  J)  - 


* 

w(0)u(t)  r 

%}  * 

2 

J 

sU ,  +  f  )/(;,  -  ~ 


— oc  — oo 


+ 


oo 

oo 

w(0)»v(t) 

J  % 

— oo  - 

-oc  —r 

(88) 


s(/,  +  ±  )s*(/1  -  f  )e>™‘'  ~  f  ^e^^dodt 


and,  by  convenient  change  of  variables, 


OC  oc 

/*  (* 


IPS(1.  J)  = 


J  J 

— oo  — oo 


+ 


u-(Q)u-(t) 

2 


m(0)u(t) 


s(^)s  (m  —  T)e/2"wV^c~'/2"<u  '^dvd r 


s  (m)s(u  +  T)ei*"t,udne  J2r{u+'^dodr  . 


(89) 


—OC  —  OC 


— oo 


Interchaneine  the  order  of  intecration,  we  obtain 


OC 

lPS(t.f)  =  y  f 


J 

— OC 


[j(/)/(r  —  t)  +  s'(/)s(i  4-  T)]u-(0)u-(r)<?~yi  '  dr  . 


completing  the  proof. 

We  can  now  easily  determine  which  properties  are  maintained  after  the  windowing 
operation,  by  direct  application  of  the  results  in  [Ref.  14),  mapping  the  characteristics 
of  the  kernel  function  to  the  properties  of  the  distribution.  The  conditions  on  d>(i>.  r) 
necessary  for  the  preservation  of  each  of  the  properties  will  also  be  given. 
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I .  Time  shift 


If  oj(t)  —  s(t  —  t0),  then 

IPSJt,J)  =  IPSs(t-t0.J). 

The  requirement  on  <T>(u,  t)  is: 

0(o,  t)  constant  for  all  t . 

This  property  will  thus  be  maintained  for  every  used  window. 


2.  Modulation 

If  co(r)  =  s(t)e:2”' 6',  then 

IPSJt,j)  =  IPSs(tJ~f0). 

The  requirement  on  0(u,  r)  is: 

0(u.  t)  constant  for  all /. 

This  property  will  thus  be  maintained  for  every  used  window. 


3.  Marginal  in  Time 

OO 

IPS(t.  f)  df=  |  s(t)  r 

%. 

—  DC 

The  requirement  on  t)  is: 

0(t\  0)  =  1  for  all  v  . 

Maintained  ill  u-(0)  =  1  . 


4.  Marginal  in  Frequency 


IPS(t ,  f)dt  =  \S[f) 


The  requirement  on  0(a.  x)  is: 

0(0,  t)  =  I  for  all  r  . 
Not  maintained  for  practical  windows. 


(91) 

(92) 


(93) 


(94) 


(95) 


(96) 


(9~) 


(98) 


30 


5.  Realness 


(99) 


lPS(t.  f)  is  real  for  all  r,  f. 

The  requirement  on  <l>(o,  x)  is: 

0(i>,  r)  =  cD*(  —  o,  —  t)  .  (100) 

Maintained  if  the  window  is  real  and  even,  or  at  least  possesses  conjugate  symmetry 
around  the  origin. 

6.  Instantaneous  Frequency 

O© 

/* 

/■ips(i ,  M 

v 

— — - - - =  fit)  (instantaneous  frequence)  (101) 

I  5(0  I 

The  requirement  on  ch(i>.  r)  is: 

0(t>.  0)  =  1  and  -r~  O(o,  r)  j.=r  =  0  for  all  u.  (102) 

CT 

Maintained  if 

•  v.-(0)  =  1  . 

•  u-'(  x  )!_..,*=  0. 

7.  Group  Delay 


J  t-IPSu.  f)dt 

— - : - =  TJfi  (group  delay)  (103) 

i  s(D  r 

The  requirement  on  <!>( i ■,  x)  is: 

<1>0>.  x)  =  1  and  x)  |, =  0  for  all  x.  (104) 

cv 

Not  maintained  for  practical  windows. 

S.  Zero  Power 

5(/0)  =  0  =*  /PS(r0.y)  =  O  (105) 

This  property  is  clearly  maintained,  due  to  the  multiplicative  nature  of  the  win¬ 
dowing  operation. 

S(fs)  =  0  =>  IPS[t.  fQ)  =  0  (106) 

This  property  will  not  be  maintained  since,  from  (  64  ),  a  convolution  along  the 
frequency  axis  will  be  performed. 
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C.  SPECTRAL  RESOLUTION 
I.  Time-invariant  Components 

Though  IPS  is  a  tool  directed  towards  the  analysis  of  signals  with  time- varying 
spectral  contents,  some  insight  into  its  inner  workings  will  be  gained  by  considering  its 
resolution  capabilities  for  time-invariant  spectral  components.  A  comparison  of  resol¬ 
ution  capabilities  lor  IPS  and  \VD  will  be  made,  and  contrasted  with  the  standard 
Fourier  Transform. 

As  seen  in  (  6  ).  IPS  is  a  coherent  average  of  two  terms,  one  of  which  uses  only 
past  data,  while  the  other  uses  only  future  information.4  When  analyzing  finite  duration 
signals,  the  "past  term"  will  have  its  maximum  resolution  capability  at  the  end  of  the 
data  segment,  since  it  is  at  this  point  that  the  most  past  data  is  available.  Similarly,  the 
"future  term"  will  resolve  better  at  the  beginning  of  the  data  segment.  This  can  be  seen 
in  Figure  6  and  Figure  7,  for  the  past  and  future  terms,  respectively. 


Figure  6.  Sinusoidal  component.  Past  term.  3-D  (a)  and  Contour  (b)  plots 


If  Tie  two  terms  are  averaged  as  in  (  6  ),  the  resulting  expression  (IPS)  will  keep  the  good 


resolution  at  end-points  that  each  component  provides,  degrading  its  resolution 


A  note  of  caution  should  be  made  here.  These  "past"  and  "future"  terms  should  not  be 
identified  with  the  two  terms  in  (  12  ),  since  each  one  of  the  terms  in  this  formula  cpans  all  the  data, 
from  —  oo  to  co. 
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capabilities  towards  the  center  of  the  time  segment.  This  effect  can  be  observed  in  the 
contour  plot  of  Figure  8. 


Figure  7.  Sinusoidal  component.  Future  term.  3-D  (a)  and  Contour  (b)  plots 


Figure  8.  IPS  for  a  sinusoidal  component.  3-D  (a)  and  Contour  (b)  plots 
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This  is  one  of  the  characteristics  of  the  un\vindovred  IPS.  It  has  better  resolution  at 
end-points  than  it  does  in  the  middle  of  the  data  segment.  We  should  point  out  that  the 
WD  does  the  opposite,  presenting  a  thinner  main-lobe  at  the  center  of  the  time  interval, 
and  loosing  resolution  capabilities  towards  the  end-points  (Figure  9). 


Figure  9.  WD  for  a  sinusoidal  component.  3-D  (a)  and  Contour  (b)  plots 


A  more  formal  discussion  can  be  made  as  follows:  Let  us  denote  by  CIF,0(t)  the 
function  of  r  whose  Fourier  Transform  is  lVD(t0,  f)  and,  similarly,  denote  by  C7,  ( r)  the 
function  of  t  whose  Fourier  Transform  is  IPS(t0,  j).  Hence, 


Clf '/„(*)  =  M'o  -  -J  H* o  +  Y  ) 

CI,b(t)  =  y  [.<r0)5*(r0  -  r)  +  s\t0)s(r0  +  r)J 


(107) 


where  s(r)  is  the  stationary  signal  of  interest.  That  is,  both  IPS  and  WD  create,  from  the 
data,  a  new  "signal",  whose  Fourier  Transform  is  the  wanted  spectrum  (Figure  10).  It 
is  thus  the  effective  duration  of  this  new  signal  that  determines  the  resolution  obtainable 
in  the  frequency  domain. 
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Figure  10.  IPS  and  YYD  :  conceptual  diagram 

Let  us  now  assume  that  s(t)  is  only  known  within  a  Finite  segment  oflcngth  L, 
and  compare  the  theoretical  resolution  of  IPS  and  YVD  at  both  the  end-points  and  the 
middle  of  the  segment. 

a.  End-points 

Since  both  end-points  are  treated  similarly,  we  will  only  consider  the  start¬ 
ing  point  of  the  known  signal,  which  we  will  denote  by  Let  us  consider  that  a 
rectangular  window  of  length  L  was  applied  to  the  signal  (  Figure  11).  It  is  easily  seen 
from  (  107  )  that 


CIV,  = 

'o 


0  T^O 
I  ■s(i’o)  I2  T  =  0 


(108) 
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and  thus,  1FZ),0(/)  (Fourier  transform  of  a  delta  function)  is  a  constant.  On  the  other 
hand,  the  windowing  of  the  data  implies  that  a  window  (  w,(/)  )  of  length  2L  is  applied 
to  C/,0( t)  (Figure  12). 


Figure  12.  Window  applied  to  C/,0(r) 
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When  we  now  take  the  Fourier  transform  of  C/,o( r)  ,  we  will  be  taking  the  transform  of 
a  signal  twice  as  long  as  the  original  data.  This  creates  a  main-lobe  twice  as  narrow  as 
the  main-lobe  that  would  be  obtained  if  we  were  to  Fourier  transform  the  data.  How¬ 
ever,  we  know  from  (  107  )  that  CllfS{r)  =  CI‘0(  —  t).  This  then  implies  that,  from  the 
obtained  segment  of  length  2L,  only  half  contains  information.  As  a  result,  at  end¬ 
points,  IPS  presents  a  main-lobe  twice  as  narrow  as  the  one  obtainable  by  Fourier 
transforming  the  data,  but  without  an  equivalent  improvement  in  effective  resolution. 
b.  Middle-point 


We  now  consider  that  the  window  is  symmetrically  placed  around  i0,  as  in 

Figure  13. 


Figure  13.  Symmetric  window  applied  to  the  signal 


Again,  it  is  readily  seen  from  (  107  )  that  windowing  the  data  with  a  window  of  length 
L  centered  at  t0  implies  the  windowing  of  C1F,0(t)  with  a  window  of  length  2L.  The  effect 
on  C/,o(t),  however,  is  to  window  it  with  a  window  of  length  L.  This  discrepancy  implies 
that,  in  the  middle  of  the  interval,  the  WD  has  twice  the  effective  resolution  of  IPS.  Let 
us  note  that,  as  discussed  before,  for  both  CfF,0(r)  and  C/,0(r),  only  half  of  the  length 
carries  information.  Hence,  their  Fourier-transform  will  have  its  main-lobe  artificially 
narrowed  by  a  factor  of  two,  without  improvement  in  effective  resolution.  In  summary, 
the  WD  has  at  the  middle  of  the  time  interval  the  same  effective  resolution  that  IPS  has 
at  end-points,  and  hence  the  same  effective  resolution  obtainable  by  Fourier- 
transforming  the  data. 


37 


This  fact  of  IPS  having  twice  as  poor  resolution  in  the  middle  of  the  time 
interval  than  at  end-points  will  be  important,  since  the  windowed  implementation 
normally  uses  the  window  centered  at  the  time  of  interest.  This  then  means  that  these 
implementations  are  using,  at  each  time,  the  worst  estimate  (in  terms  of  resolution)  that 
IPS  will  provide.  The  model  of  a  typical  windowed  implementation  of  IPS  can  be  seen 
in  Figure  14. 


Figure  14.  Typical  implementation  of  IPS 

The  overall  result  is  that,  for  these  implementations,  IPS  has,  by  a  factor  of  two,  lower 
effective  resolution  for  stationary  components  when  compared  with  the  WD. 

2.  Time-varying  Components 

One  of  the  most  noticeable  characteristics  of  IPS  is  the  fact  that  the  width  of  its 
main-lobe  depends  on  the  dynamics  of  the  signal.  A  typical  example  can  be  seen  in 
Figure  15,  where  a  quadratically  chirped  FM  signal  is  used. 
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Figure  15.  IPS  for  FM  signal  -  Quadratic  chirp.  3-D  (a)  and  Contour  (b)  plots 


This  test  signal  has  its  instantaneous  frequency  increasing  quadratically  with  time.  The 
width  of  the  main-lobe  is  directly  proportional  to  the  slew  rate.  In  Figure  16,  a  cross 
section  along  the  frequency  axis  of  1  PS  for  the  same  test  signal  is  presented.  The  past 
and  future  terms  are  also  indicated. 


Figure  16.  IPS  for  FM  signal  -  Quadratic  chirp.  Cross-section 


As  is  seen,  the  width  of  the  main-lobe  of  IPS  results  from  the  imperfect  alignment  of  the 
two  terms:  the  past  term  seems  to  lag  the  true  position  of  the  instantaneous  frequency, 
showing  a  tendency  to  be  positioned  in  locations  previously  occupied  by  the  signal,  while 
the  future  term  leads  the  true  location  of  the  instantaneous  frequency.  This  lag  and  lead 
effect  become  more  severe  for  faster  dynamics,  hence  broadening  the  main-lobe  of  IPS. 

In  an  attempt  to  explain  this  feature,  one  might  be  tempted  to  argue  that  IPS 
lacks  symmetry  in  its  definition.  Despite  being  clearly  insufficient,  as  we  will  see,  it  will 
be  instructive  to  pursue  this  type  of  reasoning  a  little  further,  since  some  additional  re¬ 
sults  will  emerge.  When  computing  C/,fl( r),  we  are,  for  each  t  (without  loss  of  generality 
assume  r  >  0  ),  coherently  averaging  two  terms.  One  results  from  the  present  and  past 
signal  history  (jr*(r)jc(r  +  t)  ),  and  the  other  results  from  present  and  future  signal  history 
(x(t)x‘(t  +  r) )  .  To  obtain  IPS,  we  will  extract  the  frequency  contents  of  the  sequence  so 
obtained.  However,  if  the  signal  has  frequency  dynamics,  each  one  of  the  terms  will 
contribute  differently  to  the  frequency  contents  of  C/,0( t)  .  For  example,  if  the  instanta¬ 
neous  frequency  of  the  signal  is  increasing  in  time,  the  contribution  of  the  present-future 
term  (.v'(t).v(/  4-  t)  )  will  have  higher  frequencies  than  the  present-past  term,  since  it  is 
centered  in  a  region  with  higher  average  instantaneous  frequency.  The  temptation  thus 
arises  to  compensate  for  the  difference  in  the  centers  of  the  two  regions.  This  effectively 
corresponds  to  forcing  the  two  terms  in  Figure  16  to  align  "correctly",  thus  eliminating 
the  broadening  of  the  main-lobe.  To  do  so,  it  suffices  to  introduce  compensating  phase 
terms  in  the  process  of  computing  IPS  (see  .  Alternate  way  of  computing  the  Wigner- 1  ille 
below).  When  these  compensating  terms  are  introduced,  the  distribution  that  results  is 
the  Wigner-Ville  distribution.  We  can  now  interpret  WD  as:  the  WD  is  the  distribution 
that  results  by  aligning  the  two  terms  of  IPS. 

•  Alternate  way  of  computing  the  Wigner-Ville 


From  the  previous  discussion,  an  alternate  way  of  computing  the  WD  can 
be  derived,  which  allows  us  to  directly  apply  the  WD  to  real  signals  sampled  at  the 
conventional  Nyquist  rate,  thus  avoiding  the  need  for  oversampling  that  would  be 
required  if  the  WD  was  to  be  computed  in  the  usual  way  [Refs.  7,  10].  The  two 
step  procedure  is: 


1.  Compute 


F( t,  »)  =  y 


[s(r)s  (/  —  T)e~jnv'  -fs  (t)s(t  +  t )e/suT]e/2'l'V/  (109) 


— oo 
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2.  Compute 


WD{t,j)  = 


F( t,  o)e  ^  "l'ht+'^dodT 


(110) 


— oo  —  oo 


As  can  easily  be  recognized,  the  need  for  a  2-D  Fourier  transform  places  the 
cost  of  this  procedure  well  above  the  one  needed  if  an  interpolation  scheme  is  used.  This 
consideration  may  disappear  in  emerging  fields  such  as  optical  signal  processing. 

All  our  previous  discussion  has  been  based  on  a  loosely  defined  lack  of  symme¬ 
try,  intimately  related  with  the  concepts  of  past  and  future.  As  pointed  out  earlier,  this 
explanation  is  clearly  insufficient.  To  illustrate  why,  a  different  definition  of  "Instanta¬ 
neous  Power  Spectrum"  is  made  in  Appendix  E,  thus  creating  a  new  distribution.  The 
important  point  to  be  made  is:  not  only  are  the  concepts  of  past  and  future  inexistent 
in  this  new  definition,  but  it  also  possesses  perfect  symmetry.  However,  the  resulting 
distribution  is  very  similar  to  IPS,  with  the  characteristic  data-dependent  width.  Since 
lack  of  symmetry  cannot  be  argued  for  this  definition,  the  explanation  for  the  broaden¬ 
ing  of  the  main-lobe  must  be  found  somewhere  else. 

In  Appendix  C,  the  issue  of  uncertainty  is  addressed.  One  important  result  to 
extract  from  there  is:  the  maximum  obtainable  frequency  resolution  when  analyzing 
signals  with  unknown  frequency  dynamics  is  yjdfjdt .  This  result  is  in  agreement  with  the 
fact,  proven  by  Rihaczek  [  Ref.  4],  that  signals  with  strong  phase  modulation  have,  at 
each  time,  their  energy  concentrated  within  a  frequency  band  Bd  of  size 


1  <TC>(f) 


dr 


(HI) 


where  O(t)  is  the  signal's  phase.  It  thus  explains  the  broadening  of  the  main-lobe  of  IPS 
in  terms  of  the  expression  of  uncertainty.  Since  the  uncertainty  region  grows  for  in¬ 
creasingly  faster  dynamics,  so  does  the  width  of  the  main-lobe.  Put  in  simpler  terms,  "the 
faster  it  moves,  the  harder  it  is  to  locate".  According  to  the  results  of  Appendix  C.  any 
algorithm  should  behave  in  this  manner,  unless  it  has,  or  assumes  (directly  or  indirectly), 
information  concerning  the  dynamics  of  the  signal.  However,  not  assuming  is  the  key 
rule  for  a  robust  algorithm,  one  whose  performance  does  not  depend  on  the  class  of 
signals  being  analyzed. 


D.  THE  CHOICE  OF  A  WINDOW 

We  know  from  (  66  )  that,  if  we  pre-window  the  signal  with  a  window  real  and 
symmetric  around  the  time  at  which  we  are  evaluating  IPS,  the  resulting  IPS  will  be: 


ips{s.wp,j)=  ipss(t,»yips,(tj-»)du 


(112) 


But 


O© 

rPSw(t ,  J)  =  -f  -  r)  +  w(t)w(i  +  T)]e~JWrdr  (113) 

— oc 


and,  since  w{t  —  t)  =  vv(f  +  t), 


IPS„{t.  J)  =  I  w{t)w{t  +  x)e  Jlrj~'dr 


=  wo(0)iro(/) , 


(114) 


where  ivc(r)  =  w(/  +  r),  that  is.  wc  is  the  used  window  shifted  to  the  origin. 

Hence,  if  a  real  and  even  sliding  window  is  used,  the  smoothing  function  in  (  112) 
is  the  Fourier  Transform  of  the  window  w0(r)  .  It  is  interesting  to  compare  this  result 
with  the  Pseudo-Wigner  Ville  (frequency  smoothed  Wigner-Ville),  where  the  smoothing 
function  is  not  the  Fourier  Transform  of  the  window  itself,  but  the  Fourier  Transform 
of  the  square  of  the  window  (in  a  rescaled  frequency  axis)  [Ref.  6].  We  thus  see  that  all 
the  knowledge  on  window  functions  for  the  Fourier  Transform  can  be  directly  applied 
to  the  choice  of  a  window  to  be  used  with  IPS. 

Though  the  smoothing  function  in  (  1 12  )  is  always  real  (we  are  considering  only  real 
and  even  windows),  it  may  be  negative.  This  can  contribute  to  the  presence  of  strongly 
negative  values  in  the  resulting  IPS,  especially  if  the  side-lobe  structure  of  the  window 
is  not  controlled.  This  thus  suggests  that  windows  with  good  side-lobe  suppression 


42 


should  be  preferred.  Spectral  resolution  is  another  important  issue  to  be  considered 
when  choosing  a  window. 

As  discussed  in  the  previous  section  (Section  III-C),  IPS  has  a  main  lobe  whose 
width  is  proportional  to  dfjdi ,  in  accordance  with  the  uncertainty  principles  of  Ap¬ 
pendix  C.  Hence,  degrading  the  apparent  frequency  resolution  of  the  unwindowed  IPS 
does  not  become  a  relevant  issue,  if  working  with  windows  whose  effective  duration  is 
well  above  the  reciprocal  of  the  uncertainty  region.  The  width  of  this  region  should  thus 
be  the  criteria  for  the  choice  of  window  size,  but  would  require  a-priori  knowledge  of  the 
dynamics  of  the  signal.  In  the  absence  of  this  knowledge,  optimality  will  not  be  achieved. 
However,  as  we  will  see.  IPS  is  a  very  forgiving  tool,  in  the  sense  that  its  performance 
is  not  seriously  affected  even  for  large  deviations  from  the  optimal  window  size. 
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IV.  EXPERIMENTAL  RESULTS 


To  illustrate  the  types  of  behavior  predicted  in  the  previous  sections,  some  test  cases 
will  be  presented.  For  all  the  cases  the  test  signals  are  128  point  sequences,  with  one 
spectrum  for  each  sample.  All  plots  are  on  a  linear  scale.  Some  considerations  relating 
to  the  Pseudo  Wigner-Ville  (PWD)  were  made:  all  test  signals  are  analytic,  thus  freeing 
the  PWD  from  the  somewhat  annoying  interferences  between  positive  and  negative  fre¬ 
quencies.  The  exception  to  this  rule  will  be  Figure  19,  since  its  sole  purpose  is  to 
illustrate  the  elTects  of  the  use  of  real  signals.  Also,  when  noise  is  present  (Figure  24), 
analytic  noise  is  used.  Unless  otherwise  stated,  all  Figures  use  a  41 -point  Hamming 
window,  and  present  the  distributions  after  smoothing  along  the  time  direction.  Hence, 
most  of  the  spectral  cross-terms  of  the  PWD  due  to  multi-component  signals  do  not 
appear  in  the  plots.  To  illustrate  the  effectiveness  of  time-smoothing,  Figure  17  was 
provided,  differing  from  Figure  21  only  in  the  fact  that  the  time-smoothing  technique 
was  not  used. 


Figure  17.  FSK.  IPS  (a)  and  PWD  (b).  No  time-smoothing. 


A.  REAL  SIGNALS 

In  Figure  18,  IPS  and  PWD  for  a  linear  analytic  FM  chirp  are  presented. 
Comparing  Figure  18  with  Figure  19,  where  the  test  signal  is  a  linearly  chirped  cosine 
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function,  we  can  appreciate  how  insensitive  IPS  is  to  the  fact  of  the  signal  being  real. 
The  artifacts  present  in  the  PWD  when  analyzing  real  signals  are,  for  all  practical 
purposes,  absent  in  IPS. 


Figure  18.  Analytic  linear  FM  chirp.  IPS  (a)  and  PWD  (b). 


Figure  19.  Linearly  chirped  cosine.  IPS  (a)  and  PWD  (b). 
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TIME  AXIS 


B.  INSTANTANEOUS  POWER 


Another  visible  ellect  in  Figure  19  is  the  "amplitude  modulation”  of  the  main-lobe 
of  IPS.  IPS  gives  us  the  distribution  in  frequency  of  the  instantaneous  power  of  the  sig¬ 
nal.  It  then  follows  that,  when  computing  the  spectrum  at  a  low  power  sample,  the 
amplitudes  of  the  frequency  components  will  have  to  be  small,  since  their  total  power 
must  add  up  to  the  low  power  of  the  sample.  This  is  a  direct  consequence  of  (  95  )  and 
is  a  requirement  that  has  be  fulfilled  by  any  positive  distribution.  It  is  not  observable  in 
Figure  18,  since  the  used  analytic  signal  has  the  same  power  at  all  samples. 

The  WD,  though  also  obeying  the  Marginal  in  Time  property,  tends  to  absorb  the 
fluctuations  of  the  instantaneous  power  in  its  cross-terms,  and,  in  order  to  do  so,  is 
forced  to  assume  strongly  negative  values. 

C.  END-POINT  RESOLUTION 


Though  the  better  resolution  at  end-points  provided  by  IPS  is  in  general  well  visible 
in  all  the  plots  that  are  presented,  Figure  20  is  probably  one  of  the  most  representatives. 


Figure  20.  FSK  (127-point  Hamming  window).  IPS  (a)  and  PWD  (b). 
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Considering  it,  we  see  that  not  only  is  the  end-point  resolution  oflPS  better,  but  also 
the  transition  times  are  shorter,  allowing  a  more  accurate  definition  of  the  time  ol  spec¬ 
tral  jumps. 

As  stated  in  Section  Ill-D,  IPS  is  a  forgiving  tool,  in  the  sense  that  we  can  deviate 
from  the  optimal  length  of  window  without  severely  degrading  its  performance.  To 
illustrate,  figure  21  is  presented,  which  differs  from  Figure  20  only  in  that  it  uses  a 
window  with  different  effective  duration. 


Figure  21.  FSK  (  51-point  Hamming  window  ).  IPS  (a)  and  PWD  (b) 

As  is  seen,  the  time  resolution  is  basically  kept  the  same,  the  only  noticeable  difference 
being  the  width  of  the  window's  main-lobe. 

D.  MULTI-COMPONENT  SIGNALS 

Performance  for  multi-component  signals  is  addressed  in  Figure  22  and  Figure  23, 
where  the  test  signal  has  three  components:  a  linearly  chirped  FM  signal,  a  quadratically 
chirped  1  M  signal,  and  a  stationary  component.  As  seen,  IPS  is  lice  of  the  cross-terms 
appearing  in  the  PWD.  However,  for  iPS,  each  one  of  the  spectral  components  has  its 
width  affected  by  the  associated  uncertainty  region  and  is,  hence,  wider  than  the  corre¬ 
sponding  one  in  the  PWD.  Again,  the  instantaneous  power  modulation  is  visible  in  IPS. 
These  effects  arc  also  apparent  in  the  provided  contour  plots. 
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Figure  22.  IPS  for  Multi-component  signal.  3-D  (a)  and  Contour  (b)  plots 


Figure  23.  PWD  for  Multi-component  signal.  3-D  (a)  and  Contour  (b)  plots 


E.  PERFORMANCE  IN  NOISE 

One  issue  still  not  addressed  is  the  performance  in  noise.  In  Appendix  I),  a  proof  is 
given  that,  for  a  signal  embedded  in  additive  bandlimited  white  Gaussian  noise,  the 
variance  of  IPS  is,  on  the  average,5  smaller  than  the  variance  of  \VD.  Also,  from  the 
same  Appendix,  we  have  that,  when  the  bandwidth  of  the  noise  approaches  infinity,  the 
improvement  in  the  variance  given  by  IPS  approaches  the  limiting  value  of  3  dB. 
Figure  24  is  provided  to  illustrate  this  theoretical  result.  The  SNR  for  this  picture  is 
5  dB,  and  it  is  seen  that  IPS  does,  indeed,  present  a  smoother  estimate  of  the  instanta¬ 
neous  spectrum. 


Figure  24.  FSK  (  SNR  =  5  dB  ).  IPS  (a)  and  PVVD  (a) 


F.  THE  CAPTURE  EFFECT 

Having  borrowed  the  term  from  communication  theory  (FM  demodulation),  our 
goal  is  to  compare  the  behaviors  of  IPS  and  PWD  when  analyzing  multi-component 
signals,  when  the  components  have  different  amplitudes.  In  Figure  25,  the  two  compo¬ 
nents  have  amplitudes  differing  by  a  factor  a  four  (12  dB).  As  is  observable,  both  IPS 
and  PWD  detect  the  weaker  component,  despite  the  fact  of  the  cross-terms  in  the  PWD 
being  already  more  energetic  than  the  weak  signal  itself.  In  Tigure  26,  the  same  two 

■*  The  variance  of  IPS  is  time  dependent.  When  we  say  that  "...the  variance  is. ..on  the  aver¬ 
age....",  we  are  referring  to  the  time  average  of  the  variance. 
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component  signal  is  presented,  now  with  a  factor  of  8  between  amplitudes  (18  dB).  It 
is  obvious  that  the  PWD  lost  the  weak  component.  The  cross-terms,  however,  still  re¬ 
main.  In  this  case,  the  position  of  the  weaker  term  could  be  inferred  from  the  cross¬ 
terms,  but  such  a  technique  would  be  impossible  to  apply  to  more  complex  signals. 


Figure  25.  Power  ratio  =  12  dB.  IPS  (a)  and  PWD  (l>) 


Figure  26.  Power  ratio  =  18  dB.  IPS  (a)  and  PWD  (b) 
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V.  PARAMETRIC  IPS 


A.  INTRODUCTION 

We  have,  from  (  12  ),  that  the  IPS  of  a  signal  is,  at  each  time,  the  Fourier  transform 
of  a  function  of  the  signal.  This  implies  that  the  wanted  information  lies  in  the  spectral 
content  of  this  function.  The  process  of  obtaining  the  IPS  can  therefore  be  reduced  to 
the  problem  of  performing  conventional  spectral  estimation  on  a  function  G{t,  r)  which 
is  derived  from  the  signal  as  in  (  II  ).  The  use  of  parametric  methods  should  thus  be 
considered,  in  an  attempt  to  improve  the  resolution  capabilities  of  IPS.  These  consider¬ 
ations  are  also  applicable  to  the  WD,  and  a  parametric  version  of  this  distribution  has 
been  considered  [Ref.  20]. 

B.  AR  MODELING  -  CONSIDERATIONS 

Though  obtained  as  a  bilinear  transformation  of  the  data,  the  function  in  (  1 1  )  can 
not  be  considered  a  true  autocorrelation  function,  since  it  is  not  constrained  in  any  way 
other  than  possessing  conjugate  symmetry  around  the  origin.  It  can  not  be  guaranteed 
to  give  rise  to  a  positive  definite  autocorrelation  matrix,  and  hence  can  not  be  used  to 
directly  solve  the  normal  equations.  We  can  however  estimate  its  autocorrelation  func¬ 
tion  and  fit  the  mode!  to  this  autocorrelation.  In  a  sense,6  we  will  be  fitting  the  model 
in  the  fourth  moment  of  the  data.  In  all  the  results  presented,  the  AR  parameters  are 
obtained  using  the  Modified  Covariance  Method.  This  method  has  been  found  to  pos¬ 
sess  good  resolution  properties,  while  alleviating  some  of  the  problems  associated  with 
the  Maximum  Entropy  Method  (MEM)  [Ref.  21].  All  plots  are  on  a  logarithmic  scale. 

C.  EXPERIMENTAL  RESULTS 

Choosing  the  order  of  the  model  is  a  sensitive  issue.  From  Figure  16.  we  know  that 
the  main  lobe  of  IPS  consists  of  two  terms,  which  do  not  align  exactly.  We  might  thus 
expect  that,  if  the  correct  order  is  exceeded,  the  extra  poles  will  tend  to  resolve  these  two 
terms,  and  each  spectral  component  will  have  two  associated  peaks  in  the  spectrum.  In 
Figure  27,  where  the  test  signal  is  again  a  single  component  quadratically  chirped  FM 
signal,  this  effect  is  well  visible.  In  the  left,  a  one-pole  model  was  used  (analytic  test 
signal).  With  the  addition  of  one  extra  pole,  we  can  see  that  the  two  terms  were  resolved 

6  Lacking  both  temporal  and  statistical  averaging.  (11)  can  not  be  considered  to  represent  the 
second  moment  of  the  data. 


5! 


(with  some  bias),  and  the  undesired  effect  is  now  present.  In  Figure  28,  a  fourth  order 
model  was  used  to  fit  a  three  component  signal. 


Figure  27.  Parametric  IPS.  First  (a)  and  second  (b)  order  models 


We  can  observe  that,  since  the  main  lobes  were  not  over-resolved,  the  signal  components 
are  well  defined,  and  can  be  located  with  greater  precision  than  what  would  have  been 
possible  with  the  non-parametric  IPS. 


Figure  28.  Three  component  signal.  Parametric  IPS.  Fourth-order  model. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  thesis,  an  alternative  tool  for  the  spectral  analysis  of  signals  with  time- 
varying  spectral  content  is  presented  (IPS).  Having  its  roots  in  the  definition  proposed 
by  Page  of  "Instantaneous  Power  Spectrum",  IPS  has  a  defining  expression  very  similar 
to  the  one  defining  the  Wigner-Ville  distribution  (WD),  Their  performance  is,  however, 
considerably  different.  The  WD  approach,  though  extensively  used  in  the  last  decades, 
suffers  from  a  number  of  short-comings.  The  foldover  frequency  of  its  discrete  fre¬ 
quency  is  located  at  n/2  ,  requiring  the  use  of  analytic  signals  or  the  use  of  an  interpo¬ 
lation  scheme  as  an  alternative  to  sampling  at  twice  the  Nyquist  rate.  Interference  when 
used  with  real  signals,  cross-terms  when  multi-component  signals  are  present,  preference 
of  linear  dynamics  and  loss  of  resolution  at  the  extremes  of  the  analysis  segment  are 
other  problems  associated  with  the  WD  [Ref.  8].  The  distribution  presented  here  is 
shown  to  have  some  advantages  over  the  WD.  These  are:  evenness  of  performance  for 
different  dynamics  of  the  signal,  direct  applicability  to  real  signals  sampled  at  the 
Nyquist  rate,  reliability  when  analyzing  multi-component  signals,  and  preservation  of 
spectral  resolution  at  end-points.  Also,  IPS  is  shown  to  behave  better  in  the  presence 
of  additive  white  Gaussian  noise  (AWGN).  As  a  disadvantage.  IPS  does  not  achieve  the 
frequency  resolution  provided  by  the  WD.  Due  to  its  robustness,  IPS  is  especially  well 
fitted  to  be  used  as  a  front-end  tool  in  non-stationary  spectral  analysis.  The  usefulness 
of  AR  modeling  as  applied  to  IPS  is  investigated.  It  is  found  to  be  a  technique  very 
sensitive  to  the  chosen  order  model,  but  providing  some  advantage  in  defining  the  lo¬ 
cation  of  each  spectral  component.  The  use  of  other  parametric  techniques  should  be 
investigated.  In  Figure  29.  we  can  see  the  FSK  signal  of  Figure  20,  this  time  processed 
with  the  following  expression: 


OO 

f  [$(0*  (*  —  T)  s  (0*s(/  +  T T\e~J2  ^'dr  . 


— oo 


(115) 


That  is,  Figure  29  is  the  imaginary  part  of  the  Rihaczek  distribution  of  the  test  signal. 
As  is  seen,  it  shows  some  natural  ability  to  detect  fast  transitions  in  the  spectrum,  and 
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seems  to  be  a  promising  area  for  future  work.  Also,  linear  combinations  between  IPS 
and  this  imaginary  part  of  the  Rihaczek  distribution  should  be  investigated. 


Figure  29.  Imaginary  part  of  Rihaczek's  Distribution  -  FSK 


Signal  synthesis  from  IPS  may  also  be  a  promising  field  of  research  for  speech  processing 
applications. 
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APPENDIX  A.  IPS  FOR  FINITE  DURATION  DISCRETE  SIGNALS 

In  this  appendix,  a  formal  derivation  of  the  expression  defining  IPS  for  discrete 
signals  (  16  )  is  presented,  as  an  alternative  to  the  direct  discretization  of  (  12  ). 

The  discrete  version  of  (3)  is,  assuming  that  the  sequence  starts  at  n  =  0: 


S~(n,  6)  =  A  s{n)e'jen 


n=0 


(116) 


Hence, 


r 


A’=-/ 


(117) 


where 


o--Y 


ck  —  ^  s(«).s  (n  —  k)  k>  0 
n=k 


(118) 


and 


ck  =  (clk)' . 


(119) 


Similarly 


v-i 

S+(n,  6)  =  A  7^  s{n)e~jen 

n=t 


(120) 
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(121) 
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S+(«,  6 )  |2  =  A7*2  >  s(n)e~J6n  ^  s\n)e/6n 
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(-V— 1— /) 


=  A  T1  Y,  tf' 


/c=— <A— 1— r) 


where 


.v-i 


n=;+A 


s(n)s  ( n  —  A) 


A  >0 


(122) 


and 

=  (cV  • 


(123) 


Now,  to  get  the  digital  equivalents  of  (2),(4)  we  must  approximate  the  derivatives  of 
(  117)  and  (  120  ),  which  we  will  do  using  first  order  differences: 


p>.  (?)  *  AT  [  |  S>,  (?)  |2  -  |  S~(t  -  1,0)  |2  ] 

-A7-V 

k=-t 


(124) 


where 


and 


tk  =  s(«)s  (n  —  k)  k  >  0 


(125) 


(126) 


In  the  same  way. 


p+(n,  0)  *  A7[  |  S+(n,  (?)  |2  +  |  S+(i  +  1,0)  |2] 


=  AT 


O'- 1 -/) 


+  -y*A 
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k=- (V— )—  ;) 


(127) 
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where 


tl  —  s  («)s(«  +  k)  k  >  0 


(128) 


and 

4  “(O'-  029) 

Assuming  the  signal  to  be  zero  outside  the  known  samples,  then  according  to  Levin's 
definition  (6): 


IPS(n,  0)  =  -±-  [p"(«f  0 )  +  p\n,  0)] 


oc 

=  [s(rt)s  (n  —  k)  +  s  (n)s{n  +  k)']e~^ek 


k=-x 


(130) 


which  is  the  wanted  result. 
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APPENDIX  B.  INSTANTANEOUS  FREQUENCY  AND  IPS 

The  proof  of  property  6  (  55  )  can  be  made  as  follows: 


J 

*— oo 


f-lPS(t,  j)df 


I  5(f)  | 


(131) 


OO  OO 


I  5(f)  I 


[s(/)s  (f  —  t)  +  s  ( l)s{[  +  t)]c  ^^'drdf 


— oo  — oo 

oo 


I  5(f)  I2 


WO*  (/  -  T)  +  5  (t)5(f  +  ?)] 


fe  j2'f'dfdr 


2  I  5(f) 


[s(f)j*(r  -  r)  +  s*(z)s(z  +  t)]  [o(t)]c/t 

n  or 


=-JL-(-^) 

2  |  5(f)  |2  v  ^  /% 

— OC 

+  — ! —  (f) 

2  I  5(0  |  W"  / 


5(/)5  (f  —  r)<5  ( x)dr  + 


5  (r)5(f  +  t)0  (t)z/t 


- - 7  5(f)5*'(f)  -  ~  —  —  S*(f)s  (f) 


4 ft  I  5(f)  I 

j 

An  |  5(f)  |2 

j 

An  |  5 (/)  |2 


An  |  5(f)  | 
[5(f)/ (f)  -  5*(f)5  (f)] 

■[-y'2Imag[5*(f)5  (/)]] . 


(132) 
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Therefore 


OO 

/» 

f-IPS{t,  J)df 

L  Imagfs'uMO] 


5(0  I2  2n  |  5(/)  |2 


(133) 


For  real  signals, this  center  of  mass  is  zero,  as  expected,  since  the  spectrum  is  conju¬ 
gate  symmetric  about  DC. 

For  analytic  signals: 

s(f)  =  m(/)^,!>  (134) 

and 

OC 

/■ipsa.  j)df 

ty _ lmag[s*(Qs  (Q] 

I  5(0  I2  2rr  |  5(/)  |2 

Imag[m('r)t’"v"?(,)[w  (t)f'i‘in  +  j<i>  (/)A?j(^)^<r'<°]l  '  * ^ 

2  nm2(t) 

Imag^nf/lnt  (0  +j<j)  (t)m*(/)] 

2-m2(t) 


and.  finally 


fU)  = 


J _ d_ 

2  7i  dt 


[0(0] 


(136) 


which  is  the  definition  of  instantaneous  frequency  for  analytic  signals.' 


T  his  definition  of  Instantaneous  Frequency  was  introduced  by  Vide  [Ref.  2j  in  1952. 
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Hence,  we  have  that,  for  analytic  signals 


* 

f-IPS{i,  f)df 

~ ~ rr — -/xo  037) 

I  5(0  I 

and  our  proofis  complete. 
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APPENDIX  C.  NON-STATIONARY  SIGNALS  AND  THE  PRINCIPLE  OF 

UNCERTAINTY 

When  analyzing  signals  with  time-varying  spectral  contents,  there  is  a  maximum 
obtainable  frequency  resolution,  determined  by  the  dynamics  of  the  signal.  This  is  in 
contrast  with  the  stationary  case,  where  the  only  constraint  on  the  obtainable  resolution 
is  the  observation  time.  The  presentation  and  proof  of  this  effect  is  the  goal  of  this 
chapter. 

A.  THE  PHYSICAL  CONCEPT 

Before  presenting  any  mathematical  results,  we  will  try  to  approach  the  subject  with 
a  brief  discussion  of  the  involved  ideas. 

Let  us  consider  our  ability  to  perform  spectral  estimation  on  a  stationary  signal. 
For  that  purpose,  we  will  assume  that  the  incoming  signal  has  been  segmented  into 
contiguous  intervals  of  At  seconds.  Our  final  goal  is  to  determine  the  spectrum  occupied 
by  the  signal  at  roughly  the  moment  of  analysis,  that  is.  within  the  last  received  segment. 

If  we  are  given  only  this  last  segment,  the  obtainable  frequency  resolution  is  con¬ 
strained  by  |  Ref.  1] 


A/> 


1 

2  At 


(1?S) 


But  if  now,  for  some  reason,  the  next  to  the  last  segment  of  data  becomes  available, 
we  will  be  able  to  improve  our  analysis  of  the  last  segment  up  to  a  resolution  of 


Af> 


1 

4  At 


(139) 


We  can  continue  the  process  and,  if  more  segments  of  data  become  available,  im¬ 
prove  our  spectral  description  of  the  signal  at  the  present  time.  The  rule  in  this  case  is: 


•  The  stationarity  of  the  signal  allows  us  to  improve  the  analysis  of  the  present  time 
by  using  data  collected  no  matter  how  distant  in  the  past.  Hence,  if  infinite  obser¬ 
vation  time  is  allowed,  infinite  spectral  resolution  can  be  achieved. 
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Let  us  now  consider  the  same  scenario,  but  with  a  non-stationary  incoming  signal, 
which,  we  will  assume,  can  be  considered  stationary  within  the  segment.  If  only  the  last 
segment  is  available,  our  frequency  resolution  is  again  constrained  by 

A/SyL-,  „40, 

and  the  only  way  to  improve  it  is  to  consider  more  segments.  But  a  qualitative  difference 
appears  at  this  point.  Segments  of  the  “distant  past”  are  rendered  useless  by  the  dy¬ 
namics  of  the  signal  unless  we  know  or  assume  its  behavior.  Considering  those  segments 
will  only  degrade  our  spectral  analysis.  The  usable  data  is  thus  restricted  to  the  most 
“recent  past",  which  implies  an  upper  bound  on  the  obtainable  resolution. 

•  For  a  non-stationary  signal  with  unknown  dynamics,  there  is  a  maximum  obtaina¬ 
ble  frequency  resolution,  which  depends  on  how  distant  the  unusable  “distant 
past”  is.  It  is  an  absolute  maximum,  and  depends  only  on  the  dynamics  of  the 
signal. 

Proving  this  statement  is  the  purpose  of  this  Appendix. 

B.  MATHEMATICAL  FORMULATION 

Let  us  briefly  discuss  Gabor's  result.  Gabor  [Ref.  1]  defined  the  "effective  duration" 
A/  and  the  "effective  frequency  width"  A/ of  a  signal  «/>( /)  by  the  following  equations 

A t  =  y'  2 jtU  -  if 

____  (.141) 

A/=  v  2 -{f-Jf  , 

where  [.]  stands  for  the  average,  and  the  moments  are  defined,  ommiting  the  argu¬ 
ment  of  i//(/).  as 


oc 

f  \p*tr\pdf 


J 


r 


J  \p  \j/dl 

— oc 


oc 

r 


^  “7 ~  4tdf 

at 


\p  tydi 


—  yo 


(142) 
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resulting  in 


ArAf  >  —  . 


(143) 


This  only  states  that  any  signal  whose  effective  duration  is  At  (  see  Figure  30  ),  will 
have  an  elTective  bandwidth  Af  (  see  Figure  31  )  greater  than  ;  the  issue  of  uncer¬ 
tainty  is  not  addressed  in  this  result. 


Figure  31.  Effective  bandwidth  of  s(t) 


In  fact,  its  relevance  as  an  uncertainty  principle  only  comes  into  play  when  we  realize 
how  classical  Spectral  Estimation  is  achieved. 

Classical  methods  are  based  on  the  “measure  of  similarity”  between  the  data  and 
more  elementary  functions  of  time  which  have  strong  affinities  with  the  physical  concept 
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of  frequency.  The  best  example  is  the  Fourier  Transform,  which  conceptually  is  a  bank 
of  correlators,  trying  to  determine  the  “measure  of  similarity"  between  the  data  and 
sinusoidal  functions,  the  true  representation  of  our  concept  of  frequency.  Other  typical 
examples  are  in  the  works  of  Gabor  [Ref.  1]  and  Priestley  [Ref.  22],  using  elementary 
functions  on  the  i—f  plane,  where  these  functions  have  what  Priestley  calls  an 
“oscillatory’  form”  needed  to  establish  the  correspondence  to  our  notion  of  frequency. 

Now,  if  the  data  is  of  length  T,  the  elementary’  functions  can  have  at  best  the  same 

duration,  which  implies  that  their  frequency  width  is  greater  than  — r.  Since  the  ob- 

1 

tainable  frequency  resolution  is  at  best  of  the  size  of  the  frequency  width  of  the  ele¬ 
mentary  functions,  the  principle  of  uncertainty  in  Spectral  Estimation  follows. 

Gabor's  result  is  only  concerned  with  stationary’  signals,  or  at  least  assumed  to  be 
stationary,  since  A/  is  the  effective  width  of  the  Fourier  Transform  of  the  signal.  How¬ 
ever,  some  implications  can  also  be  derived  for  the  non-stationary  case. 

Let  us  assume  we  have  an  arbitrary  stationary  signal  \j/(t ),  which  we  will  consider  to 
be  analytic  to  simplify  the  formulation,  with  effective  duration  At.  effective  bandwidth 
A f  and  constant  instantaneous  frequency.  Formally, 

-r—  -j-  [  arg[t//(r)]  ]  =  constant  (144) 

ctl 

Let  us  now  force  the  signal  to  linearly  chirp  in  frequency,  by  multiplying  it  with  the 
phasor  c  T*'2  and  defining  the  modified  signal  ^w(/)  by 

^\f(t)  =  'l>(t)e''h  (145) 


where  A  is  a  positive  real  number.  Due  to  the  exponential,  i j/M(t)  is  now  a  non-stationary 
signal  chirping  linearly  in  frequency  with  a  slew  rate  of 


The  following  results  can  be  established  for  ipM(t)- 


(146) 


•  Effective  time  duration 

The  effective  time  duration  will  be 

A t»  =  v  2n(r  w  -  iMf  (147) 
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where 


'I'm'Pm  1  dt 


— oc 


\p  \p  t  dt 


l  M  oo 


*/ 

— OC 


=  t 


[ft  \ftdt 


— o© 


(148) 


Hence, 


< I'm'J'm  1  dt 


\p  \p  t2  dt 


2  \  —  oo 


(>h)  - 


-('T 


\ft  \ ftdt 


— oo 


(A/,,)2  =  2n  [  (/2)  -  (r)2]  =  (At)2 . 


•  Effective  bandwidth 

The  elTective  bandwidth  will  be 


4 f\i  —  \  2n(J.\t -JxiY  —  2^  [  (fyjf  —  if,)1  ] 


where 


(149) 


(150) 


(151) 


'Pu—'fvfdl 


t/'.v/i \>xtdl 


— OO 


Using  the  identitv 


d'"'<  .  A-  |VWf/*‘-’]  ,  tif  «**■’  +j2nk,^k'‘ . 


dt  dt 


(152) 


(153) 
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we  have  that 


1> 


M 


efk^  +  jipinkt^ 
dt  J 


di 


fst  ~ 


1 


2  nj 


ip  pdt 


— OO 


OO  OO 

r  r 

d+  V 
*  —dt 


1 


ip  \pktdt 


2  nj 


■  + 


— OO 


ip  ipdt 


(154) 


\p  \jjcJ[ 


— OO 


r  (* 

,♦  dp  . 


Ip  Iptdt 


2 */ 


■  +  k- 


ip  pdt 


ip  ipdt 


■-f+kt. 


Hence,  all  we  are  missing  to  completely  determine  AfM  is  the  expression  for  (/)',) 
where  [Ref.  1] 


d'p  \f  dipM 
dt  dt 


dt 


if  si)  = 


(2tt)2 


(155) 


'I'M'l'Sfdl 


— OO 


Noting  that 

dipXf  d>pX}  dip  dip  .dip  dip  ■>  ■>  . 

dT-*2**  — +  /2 ***-£- +  (2**)VW  ,  (156) 
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we  have 


foj)~ 


J 

— oo 


dp  dp 
di  dl 


-jlnktp  ——  + j2nktp  — ^ — h  (2nk)2t2pp 


dl 


>M)  -  n  \2 
(2  n) 


p  \pdi 


— OO 


OO 


(157) 


p  rpdt 


-(/V* 


2  — oo 


.  dp 


jk  t 


•  dip 


OC'  1 

OO 

f 

27T 

OC 

r 

p  pdt 

% 

(A 

1  p  pdt 

Considering  now  the  following  identities 


we  have 


*'i!r  m  ]+j^''*[ar8W] 


(158) 


pip  —  [  arg(0)  ]<* 


C/a/)  =  (n  +  +  IF 


A-  * 


\Jj  \jjdi 


— oo 


and,  using  (  144  )  we  get 


lti)  =  {f2)  +  ky)  +  2klf. 


(159) 


(160) 
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Since  [Ref.  2] 


'I'f'l’dt 


f  = 


In 


\p  \jjdt 


*> 

— OO 


we  have,  by  (  144  )  that 


OO 

f  dt 


f=fi—^ - =/■ 

f  \jj  \Jjdt 


— OO 


From  where. 

(^)  =  (^  +  r1(?)  +  2^y. 

Hence,  from  (  151  )  and  (  154  ). 

4A/2  =  2r[(T)  +  *V)  +  2kif  ]  -  2 n[(j)2  +  k\t)2  +  2klf  ] 

=  2^[(^-(/j2]  +  A-2[(?j-(7)2] 

=  (A/)2  +  k7(Atf  . 


(161) 


(162) 


(163) 


(164) 


Therefore,  using  (  150  )  we  can  write  the  effective  time  and  frequency  widths  of  the 
modified  signal  t//w(/)  in  terms  of  the  effective  widths  of  the  original  stationary  signal 
i //(/)  as 


(A/v,)2  =  (A/)2  +  A2(Ar)2 
(A/a/)2  =  (Ar)2 . 

Using  (  143  )  and  (  165  ),  we  get  the  result 


(165) 


f 
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(A/w-A/.v)2  =  (At)2  [  (A/)2  +  A2(A/)2  ] 
=  (ArA/)2  +  A2(A/)4 


from  where 


(A/u)2^— ! -r  +  *W 


Hence,  if  we  minimize  the  right  hand-side  of  (  167  ),  we  will  find  a  lower  bound  on  the 
effective  bandwidth  of  the  chirped  signal.  Proceeding  with  the  minimization,  we  have 


e 

cAt 


4-  k2(At)2 


_  •> 

=  0  => - - 

4(At)' 


•+2A‘(A/)  =  0 


>k2  Ar  = — - 


,  .  .4  1 

>(Af)  = - : 


/  1 

►Ar  =  ■'  — 
\  2A 


Hence 


(a/ui2>4+4 


A f.u  ^  %  k 


and  we  can  read  the  obtained  results  as  statine: 


Lemma  1.  A  linear  chirp  will  always  occupy  an  effective  bandwidth  greater  than 
N  d/jdi ,  independently  of  its  effective  duration. 
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Lemma  2.  A  linear  chirp  will  occupy  the  minimum  effective  bandwidth  of  N  dfjdi  , 
onlv  if  its  effective  duration  is 


These  results  will  have  obvious  implications  when  classical  spectral  analysis  of  the 
chirped  signal  is  attempted;  no  matter  what  the  duration  of  the  elementary  functions  is, 
the  signal  will  always  appear  to  occupy  more  than  Jdfjdt ,  and  this  will  thus  be  the  best 
obtainable  resolution.  This  completes  our  proof. 


Lemma  3. The  maximum  obtainable  frequency  resolution  when  analyzing  signals 
with  unknown  frequency  dynamics  is  N  dfjdi  . 


Though  the  previous  discussion  was  restricted  to  the  case  of  a  linear  chirp,  the 
underlying  ideas  remain  valid  for  laws  of  variation  other  than  linear;  tire  treated  case  can 
be  considered  as  a  first  order  approximation  for  arbitrary  non-stationary  signals. 
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APPENDIX  D.  PERFORMANCE  IN  NOISE 

As  we  will  see  later  in  this  Appendix,  an  absolute  evaluation  of  the  performance  of 
IPS  in  noise  is  not  possible,  because  the  second  moment  is  data  dependent.  We  will  thus, 
instead,  compare  the  noise  performance  of  IPS  with  the  noise  performance  of  the 
Wigner-Ville  distribution.  Let  us  consider  the  inverse  Fourier  transforms  of  IPS  and 
WD. 


CIs(t,  t)  =  ~  [s(r)s‘(/  -  r)  +  s*(t)s(r  +  t)] 
versus 

civs(,'T)  =  S(t+  tM'-t), 


(172) 


Given  a  signal  z{i)  that  is  corrupted  by  WAGN,  where  the  noise  is  analytic  and 
bandlimited  to  the  maximum  frequency  (IF)  of  the  signal,  that  is: 


>•(/)  =  r(/)  +  «(/) , 


(173) 


where 


S„(w) 


P'o 

to 


0<  u<  IF 
otherwise 


■> 


(1741 


(175) 


A.  MEAN 

1.  WD 


CU  v(i,  t)  =  v(t  +  —  )v  (/  —  —  )  =  £  z{i  +  ~  )  +n(t  +  —  (t  — )  4-  n  (/  — y  ) 
—  z{t  +  —  )r  (t  — )  +  z(t  +  )n  {[  — ) 

a-  i-  >- 

+  n(l  +  —  )n  (/  — —  )  4-  n(i  +  —  )r  {t  —  —  ) . 


(176) 


Hence. 


£{Cirv(t,  t)}  =  Cf1’z(r.  t)  -hz(r  +  f  )E{n\t -y)} 

+  z'0~Y  )E{n(t  +  y)«*(/~y)} 

£{arv^r)}  =  ar2(i,t)  +  /?n(T). 


2.  IPS 


Hence. 


CIA,  t)  =  y  [  v(/)v*(r  -  x)  +  v  (r)v(r  +  t)  ] 

=  yOM  +»(f)]-tV-T)  +  nV-r)] 

+  i  +'J*(/)]'[Z(/  -  T)  +  «('  -  t)]  . 


E{CIv{t ,  t)  }  =  y  r(/)r*(/  -  t)  +  y  £{«(/)«*(/  -  r)J 
+  4*  2*(0z(f  +  t)  4-  4*  £{«*(/)«(/  +  r)| 
=  CU',  t)  +  y  Rn(r)  +  y  Rn(r) . 


E{CIv.(t ,  t) }  =  CI.(i,  t)  +  . 


Therefore.  IPS  and  \VD  have,  as  far  as  their  first  moment  is  concerned,  identical  be¬ 
havior  in  the  presence  of  noise. 

B.  VARIANCE 
1.  VVD 


Civ,u,  T)  I2  =  v(r  +  4-  )v*(r  -  4-  )/(i  +  4-  )v(t  -  4- ) 

—  [’('  +  y )  +  «('  +  4'  )XZ  ('  —  4" ) +  "  ('  ~ 

■[z  (r  4-  4"  )  +  «  ([  +  4"  z  _  “7  )  +  n(t  -  4~  )] 
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(183) 


!  t)  I2  =  [z(/  +  y  )z*(f  —  f )  +  z(/  +  f  )n\t  -  y ) 

+  «('  +  y  )**(*  -  y )  +  +  y  )«*(/  -  -7 )] 

•[z*(/  +  y  )z(f  -  y  )  +  +  y  )«(■ t  -  y  ) 

+  n*(f  +  -i-)z(i“y)  +  «V  +  y)«(r-y)]. 

Using  now  the  expectation  operator,  we  have: 


E{\  Cn\(t,  t)  |2}  =  |  z(r  +  y )  |2  j  z(r  -  y )  |2 

+  z(t  +  y  )z  (/  -  y  )£{«  (t  +  y  )| 

+  |z(t  +  ^-)|2£{|»(/-y)|2} 

+  z(r  +  y  )z(r  -  -±-  )£{/»‘(t  +  )«*(r  -y)} 

+  z  (r  -  "7  )z  (/  +  y  )£-j»(t  +  y  )«('  -  —  )j 
+  |z(t-f  )|2£{|a«(/  +  -|-)|2} 

4-  E^n(t  +  -E-  )n  (t  -  )n  {i  +  ~  )n(t  -  ~  )  j 
+  z‘(r  +  -7- ): u  -  y  )£'|«(/  +  -E-  )n\t  -  -7- )  j 


(184) 


And,  since8 


£-[«*('  +  )"*(  ~jT  )j  =  £•["('  +  -7-  )'J(i  -  y )}  =  0  . 


(1S5) 


£{  |  r)  |2  }  =  |  z(t  +  ±  )  |2  i  z(f  -  -i- )  I2  +  2(1  +  \  )z\t  ~  4  )Rn( r) 


•V„M 


+  !  2(/  +  y )  I2  ~yy  +  +  y  M*  -  y  )*„(r) 


+  |  z(t  -  t. 


,2  Ao» 


+ 1  «„(t)  r  + 


'  An  If 


See  Appendix  F. 
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where  we  made  use  of  the  fact  that,  for  zero  mean  Gaussian  random  variabies 
[Ref.  23], 


•  E{: *,*2*3}  =  0 

•  £'{x1x2.v3.va)  =  E{xlx2}-  +  £{x,.xy  Efjcprj}  . 


Hence, 


£{|Cfrv(r,r)  |2}  =  |z(r  +  ^-)z*(r-^-)|2  +  ^-[|z(r  +  ^-)|2  +  |z(/-^-)| 

+  ^(T)^  (t  +  y  )r(f  —  y )  +  Rn(r)zU  +  y  )z  (f  -  y )  ( 1 87) 


+  I  ^n(T)  I  + 


•V*' 


Now, 


Var{Cirv(r,  r)}  =  E{  |  ClVy(i,  r)  |2j  -  E{CU'v(r,  r))E{ClVv(t,  r )}* 

=  E{  |  CU\U,  T)  I2}  -  [CH'-U.  T)  +  ^T)]  [Cfr;(r,  r)  +  /?>)] 
-  E{  I  CU\.(J,  r)  I2}  -  I  ;(/  +  f  )z'd  ~f)l 2  -  ]  EJ  T)  I2 
—  Cll\(r,  t)/?*(t)  —  Cll'*(r.  x)Rk(t)  . 


From  where,  using  (  1S7  )  and  (  172  ), 


Var{arv(t,  r)}  =  -  r  |  Z{1  +  -T  )  |2  +  |  z(i  -  -f 

A.VC  L  A.  A- 


+ 


•Vq»- 


(189) 


2.  IPS 

I  C/V(f,  T)  I2  =  y  I>(0*'  U  ~  t)  +  v*(r)v(/  +  t)][v*(/)v(/  -  t)  +  v(/)v*(r  4-  t)] 

=  y  M'K  (/)v*(f  -  t)v(/  -  t)]  +  y-  [v*(/)v(t  +  t)v’(r)v(r  -  r)]  (190) 

+  y  0('K*('  -  T)i-(r)v*(r  +  r)]  +  y-  [v‘(t)v(r  +  t)v(r)i*(r  +  r)]  . 

For  simplicity,  we  will  treat  each  term  in  (  190  )  separately,  and  use  the  fact  that 
the  second  and  third  terms  are  complex  conjugates  of  each  other. 


74 


•  Firs!  term 


£  j  -J  [v( /)v*(/)v*U  -  t)v(/  -  t)]  |  = 

=  ~  £{[r(/)  4-  «(/)](>  (!)  4-  «*(/)][’*(/  -  r)  +  n\i  -  r)][r(/  -  t)  +  n(t  -  r)]  ] 

=  ~J  (0  +  z(')»*(0  +  z'UWi)  +  «(')«*(')] 

•[r  (!  -  r)z(t  -  r)  +  r  (/  -  r)n(t  -  r)  (191) 

+  :(t  -  t )n  (i  -  t)  +  n  (t  -  r )«(/  -  r)]} 

=  ±  {:{!):' (i)z'(l  -  r)z(i  -  r)  +  |  z(i)  \2E{  |  n(i  -  t)  l2} 

-r  r(/)r  (/  -  r)F.[n  (, t)n(t  -  r)j  4-  z\t)z(t  -  r )E{n(t)n  (t  ~  r)} 

+  j  z(i  -  r)  rzr(  |  n{i)  I2}  4-  E\n(t)n  (t)n(t  -  x)n  (t  -  r)|]  . 

Simplifying. 

£'{  -J  [v( »•.*(/);•*(/  -  t)i(/  -  t)]  j  =  -j  1  r(f)--*(/  -  t)  |2 
+  -J  [  I  z(!)  i2  +  |  z(i  -  r)  ]2]  +  i-  z(i)z'(i  -  t)RJt)  (192) 

-  -j  z’(!)z( !  -  r,R„(  t)  +  j-  (  j  +“|  Rn(  r)  . 


•  Second  term 

F\  "J  +  T)v*U)r(f  -  rl]  j  == 

|  .  . 

=  —  U)  -  n  ! ! )][.'( /  +  r)  +  nit  +  r)][r  (/)  4-  n  (/)][r(f  -  r)  +  n(t  -  t)]  } 


(/!-(/  -  r)  z\nn(i  4-  t)  +  n\t)z(i  +  r)  +  n'u)n[!  4-  t)] 


[:  i  -  -)  -  r  (/)//(/  -  r)  4-  n  (i)zU  —  t)  4-  n  (l)n[l  -  r)] 

=  “  l/;:  (/):!/- t);U  -  t) -r  ;  (/|r(/-  z)E[n*U)n(!  -  t|] 

•  [!::'!  -  rl E'jiU  ~  z\n  (/)■;  (t\zu  -f-  r )£-[//  (i)n(t  -  *j] 

—  r  in.-: !  -  7)1.  n  Uini!  4-  r)  J-  4-  F.\n  U)n  (l)nli  4-  t  )nU  —  r)  J  j . 

r  |  t  )  - 

A  j  —  [:•  a, -.a  ~  inn/ -  r)]  \  -  ~  z'n)z  *(/).-(/  -  r)r(/  -  r) 

-  4-  4-  -p:*(fi:(  /  -  (194) 

“4  *4 


-r  •-  [/!.-/  f  r r J  -  —  :  (nr(7  -  rl/\,.(rl  -  /j  fr 
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•  Third  term  Since  this  term  is  the  complex  conjugate  of  the  second  term,  we 
have: 

£  j  -J  l>WlV  -  +  t)]  j  =  ~  z(t)z{t)z\t  +  r )z\t  -  t) 

+  -j  z(t)z\t  +  z)Rn(z)  +  -±-  z(t)z\,  -  z)Rn(z)  (195) 

+  -j  z{t)z{t  +  z)Rn(z)  +  ±  z(/)z*(r  -  t)^(t)  +  -i-  \  RJz)  |2 . 

•  Fourth  term 

£{^[v‘(0v(/  +  t)v(/)v*(/  +  t)]  j  = 

=  ~  £{[2*(/)  +  «*(/)]  0('  +  r)  +  n(t  +  T)][2(/)  +  «(/)][2*(f  +  T)  +  n\t  4-  t)]  } 

=  —  £{[z*(/)z(r  +  t)  4-  z\t)n(t  +  r)  +  n\l)z{l  4-  r)  4-  n\t)n(t  +  t)] 

■[z(f)r  (i  +  t)  +  z(l)n  (t  4-  t)  -f-  u(t):  (t  4-  r)  4-  n(t)n  ( t  4-  r)] 

=  —  [z  (i)z(r  4-  t)z(/)z  (t  +  t)  +  z  (/)-(/  4-  (/  +  f)j 

+  I  ={t)  \2E{  |  n(t  4-  t)  I2}  +  |  :(t  +  t)  |2£{  |  n(t)  |2} 

+  z(l)z  (t  +  t)£[/;  (/)«(/  +  t)]  4-  £'{/?  (. t)n(l)ii(t  +  r)n  (t  4-  r)]]  . 

-J  [v*(/)Wy  +  t)v(/)v*(/  4-  T)]j  =  -j  |  z*(r)-U  4-  r)  |* 

4 — :  ^ - i  !  r(  t)  4-|z(/  +  t)!‘’I4-  —  z  (t)z(t  4-  T  )R,A  t)  (19.) 

•4  s,  71  J  -4 

1  *  1  /  -V':  11  \’  1  i 

+  ~  -(f)-  (l  +  t)£„(t)  +  —  (  —  )  +  —  |  R,A z)  . 

Adding  the  four  terms  (  192  ),  (  194  ).  (  195  )  and  {  197  ).  we  thus  obtain  the 
expression. for  the  second  moment  of  IPS.  Now, 

Var{C7,(r.  r)j  =  £{  |  CIv(t.  t)  |2}  -  E{CIv(t.  r )'fE{CI,U.  t)l‘  .  (19M 

where 

£{c/,(t.  r)J-£{ C/Vir.  r)}*  =  [C£(/.  t)  +  Rn(z)yiCL(r,  r)  +  <(r)] 

=  CI.it,  r )CI.(t,  t)  4-  CI.it.  r)R„{r)  (199) 

4-  CI.(t,  z )Rn{z)  +  |  Rn(z)  !‘  . 


-6 


Substituting  in  (  19S  ),  we  get  that 


Var{CIJ(t,  r)}  =  E{  |  Cljt,  r)  |2}  -  E{CIVU,  r )}-E{CIv(t,  r)}* 

i  ,v0n 


4  2rr 


[  |  z(t)  |2  +  |  z(t  -  r)  |2  +  |  z(i)  |2  +  |  z(t  +  t)  |2  ] 


+  -jz(/)2*(/  -  t)R*(t)  +  -jz\t):{t  + 

+  -J  z\l)z{t  -  t)R„{t)  +  -J-  z(t)z\l  +  t)R„{ t) 

1  2  1/  AV*' 

+  yl  w i  +y(  “fr 


(200) 


We  can  now  compare  the  variances  of  CI,(t,r)  and  CH  ,(t,  t)  .  First  of  all  we 
realize  that,  in  both  cases,  the  variance  is  a  function  of  /  and  t.  If  t  =  0  ,  we  have: 


:Y0>F 


Var{C/v(/.  r)}  =  |  _(t)  f  + 


An  IF 


(201) 


a;  if 

Var{C/Fv(t.  r)}  =  — = —  ! 


(202) 


That  is. 


Var(C/,(t.  r)}  =  Yar{CIF,(/.  r)}  (r  =  0) 


(2ii3) 


For  r  0,  we  have  from  (  200  )  that 


Var{C7k.(/.  r j }  <  -j  -  [  1  r(f)  |2  +  I  z(t  -  r)  |2  +  !  z(t)  |2  +  1  z(t  +  t)  |2  ] 
+  y  I  RJr)  I  [  I  z(i):'(i  -  t)  |  +  |  z'(i)z(t  +  r)  | 

+  |  z*(t)z(t  -  r)  |  +  |  z{i)z\l  +  t)  |  ] 


(201) 


77 


Var{C7v(/,  t)}  <y~^[  |  z(t)  |2  +  j  z(i  -  r)  |2  +  |  z(i)  \ 2  +  |  z(r  +  r)  |2 

+  |  z(/)  ||  z{t  +  T)  I  +  I  z(t)  |  |  z(t  -  r)  I  +  I  :(t)  |  |  z{t  -  t)  I  (205) 


+  I  z(l)  I  I  +  T)  I  ]  +  ^  — 


'  An» 


■  ,  2 


2n 


The  time  dependency  of  the  expression  stops  us  from  a  direct  comparison  with 
t).  However,  since  the  signal  is  arbitrary,  we  have  no  reason  to  assume  any  or¬ 
dering  of  its  magnitude  at  different  points  in  time.  We  can  thus  average  over  time,  and 
get  the  following  result: 

<  Yar{C7v(/,  r)}  >  av  <  [4  <  |  z(t)  |2  >  av.  +  <  I  z(r)  |  |  z(t  +  x)  \  >  av 

+  <  I  z(t)  |  |  r(r  —  t)  |  >  av  +  <  |  z(t)  1  |  z(r  —  t)  |  >  flv 

+  <  |  z(t)  |  |  ZJ  +  T)  |  >  av]  +  {  \  (2U6) 


in 


< 


•Yr,tr 


<  u-wr>av  + 


X0u- 


av  '  2  r 


Using  (  1 S9  ). 


and  finally,  we  have 

•  T  =  0 


•  r  #  0 


<  Var {C7v(r.  x)}  >  av  <  <  Yar{Cirv(/,  r)}  > , 


(20- 


Var{C7v(/.x)}  =  Yar{CIFv(/,x)} 


<  Yar{C7„(/,  r)}  >  av  <  <  |  Var{CM'v(/.  x)}  |  >, 


COS) 


(209) 


For  the  case  x  #  0,  H>1,  we  can  make  the  following  approximations,  using 
(  189  )  and  (  200  ): 


1  /  -V01F  \  2 
<Var{C/l(rIr)}>flV*  4  '  ' 


<  Yar r)}>, 


•v„n-  \2 


CIO) 


cm 
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That  is,  when  the  bandwidth  of  the  noise  approaches  infinity,  the  difference  in 
noise  performance  between  IPS  and  WD  approaches  the  limiting  value  of  3  dB. 


APPENDIX  E.  A  DIFFERENT  DEFINITION  OF  "INSTANTANEOUS 

POWER  SPECTRUM" 

In  accordance  with  the  discussion  in  Section  III-C,  the  following  definition  of  In¬ 
stantaneous  Power  Spectrum  can  be  made. 


•  Definition  -  Instantaneous  Power  Spectrum  at  the  time  of  a  given  sample  is  the 
frequency  content  of  the  stationary  spectrum  implied  by  the  existence  of  that 
sample. 

That  is:  reading  the  Fourier  transform  as  a  change  of  basis  in  an  N-dimensional 
space,  the  Instantaneous  Power  Spectrum  will  be  the  difference  in  coordinates  (referred 
to  the  phasor  axis)  between  the  location  of  the  signal  in  the  \-dimensional  space  and  its 
location  in  the  (X-I)-dimensional  space  that  results  from  suppressing  the  dimension 
corresponding  to  the  sample  of  interest. 

Since  in  terms  of  the  internal  product,  the  absence  of  any  dimension  can  be  repres¬ 
ented  by  a  zero  component  of  any  of  the  vectors  along  that  dimension,  this  definition 
leads  to  a  three-step  procedure: 

1.  Compute 

I  A'M  |2 ,  (212) 

where 

or 

A’M  =  ^  x(n)e-JMn  .  (213) 

/?=— oo 


2.  Define  a  new  sequence,  obtained  by  zeroing  x(n)  at  the  sample  of  interest,  that  is: 

(214) 


■V»  -  {0t(n) 


n=h  r 
n  —  r 


where  r  is  the  time  at  which  the  Instantaneous  Power  Spectrum  is  to  be  determined, 
and  obtain 


A'-M 


(215) 


3.  Obtain  the  Instantaneous  Power  Spectrum  as: 

2 

Instantaneous  Power  Spectrum  =  |  X{c>)  |  —  |  ,\',(w)  j  .  (216) 


SO 


The  resulting  distribution  will  be  (  using  extended  sequences,  as  in  Appendix  A  ) 


|  A\co)  !2  =  A'(to)X\a>) 


=  J  x{n)e  Jcun- 

n=— oo 


and,  after  some  manipulations, 


|  X{a>)  |2  =  Y,Cke~J*k  . 

— oo 


where 


OO 

ck  =  (n  —  k)  k>0 

n=k 


and 


c-k  ~  ck 


(217) 


(218) 


(219) 


(220) 


Correspondingly. 


where 


and 


|  |2  =  Ybke-ja,k  , 


bk  =  Yx2(n)x2(n  —  k)  k>  0 

n—k 


(221) 


(223) 


SI 


Hence,  we  have  that 


oo  oo 

A'M  |2  -  |  X2(co)  I2  =  £  cke~j(uk  -  £  bke~jwk 

i fe=—oo  k——or> 

=  £  Ake-jMk-\x(r)\2, 


where 


Ak  =  [jr(r).v*(r  —  k)  +  x(r  +  k)x  (r)]  k  >  0 


A_k  =  Ak.  (226) 

that  is,  this  definition  leads  almost  exactly  (up  to  a  constant,  for  fixed  time)  to  IPS. 
What  is  important  to  realize  is  that  this  definition  has  no  concept  of  past  or  future,  and 
that  all  points  are  equally  treated,  no  matter  what  point  they  occupy  in  the  sequence. 
This  result  seems  to  put  an  end  to  our  hope  of  finding  assymetries  in  the  IPS  concept. 

As  a  final  note,  a  more  formal  proof  of  the  relation  between  this  distribution  and 
IPS  can  be  made  as  follows:  Define 


!l 

£ 

n  #  r 

n  =  r 

(22“) 

(0 

n  r 

-  |.v(n) 

n  =  r 

(22S) 

Hence, 


jf(«)  =  x.(n)  +  d(n) . 


A'M  =  F{x(n)} 
D(co)  =  F{ci{n)) 

A'M)  =  fW»)}  • 


82 


Then, 


X(o)  =  A\(gj)  +  D(aj) 


(231) 


and 

|  X((o)  |2  -  |  X7((o)  | 2  =  (A\M  +  D{co))(XI(cj)  +  />*(«))  -  Xz{w)X]{u) 
=  A,MZ)*(a))  +  A'*(co)D(w)  +  D((o)D{io) . 

But  we  know  that  D(co)  =  x(r)e~ Therefore, 

|  A'(co)  |2  -  |  X7(co)  |2  =  2Re{A'r(<y)/(/V‘°'}  +  D(co)x  (r)^ 

=  2Re{  [A'»  4-  Z>(co)]  x*(r)^r}  -  |  x(r)  |2 
=  2Re{  A\w).r*(r)^W}  -  |  x(r)  |2 
=  2- A t  IPS~  |  .r(r)  | 2 . 


(232) 


(233) 
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APPENDIX  F.  ON  THE  SECOND  MOMENT  OF  ANALYTIC  SIGNALS 

The  purpose  of  this  appendix  is  to  proof  the  following  theorem: 

Theorem:  For  an  arbitrary  analytic  random  signal  z{t),  defined  as 

z{t)  =  jr(r)  +jx{t) ,  (234) 

where  x(t)  is  the  Hilbert  transform  of  x(t),  the  following  identity  is  satisfied: 

£{2(«)z(s)}  =  0  (235) 


Proof: 

E{z(n):(s)}  =  £{  [*(«)  + jx(n)f[x(s)  +jx(s)1  } 

(236) 

=  £{.y(»).v(5))  +7'£{.Y(/7).r(5)}  +./E{.r(«)*(s)}  -  £{.y(w).v(5)}  . 

But,  from  the  properties  of  the  Hilbert  transform,  we  have  that 

Rri(?)  =  -Rix(r) 

(237) 

Rxx(~)  =  Rxx(r) 

Hence. 

£{z(/7)r(s)}  =  Rxx(n  —  s)  +  jRxx(n  —  s)  4-  jRxx(n  —  s)  —  Rxx(n  —  s) 

(23S) 

=  0. 

Our  proof  is  thus  complete. 
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